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A brief review of the evolutionary progress in computational aerothermo dynamics is 
presented. The current status of computational aerothermodynamics is then discussed, 
with emphasis on its capabilities and limitations for contributions to the design process of 
hypersonic vehicles. Some topics to be highlighted include: (1) aerodynamic coefficient 
predictions with emphasis on high temperature gas effects; (2) surface heating and tem- 
perature predictions for thermal protection system (TPS) design in a high temperature, 
thermochemical nonequilibrium environment; (3) methods for extracting and extending 
computational fluid dynamic (CFD) solutions for efficient utilization by all members of 
a multidisciplinary design team; (4) physical models; (5) validation process and error 
estimation; and (6) gridding and solution generation strategies. Recent experiences in 
the design of X-33 will be featured. Computational aerothermodynamic contributions 
to Mars Pathfinder, METEOR, and Stardust (Comet Sample return) will also provide 
context for this discussion. Some of the barriers that currently limit computational 
aerothermodynamics to a predominantly reactive mode in the design process will also be 
discussed, with the goal of providing focus for future research. 


Introduction 

A erothermodynamics couples the disci- 
plines of aerodynamics and thermodynamics. It 
most often deals with problems in hypersonic flight in 
which high temperature gas effects strongly influence 
the fluid forces (pressure, skin friction), energy flux 
(convective and radiative heating), and mass flux (ab- 
lation) on a vehicle. Hypersonic flows are usually char- 
acterized by the presence of strong shocks and equi- 
librium or non-equilibrium gas chemistry. Accurate 
prediction of these effects is critical to the design of 
any vehicle which flies at hypersonic velocities. Fluid 
forces are integrated over the complete configuration 
to define the aerodynamic forces (lift, drag, pitching 
moment, control surface effectiveness). Peak temper- 
atures, peak heat transfer and heating load (heating 
rate integrated over time) are mapped over the vehi- 
cle surface to design the thermal protection system. 
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Pressure distributions are required for assessment of 
structural loads and venting environments. 

Computational aerothermodynamics uses numerical 
solutions of the governing equations for continuity, mo- 
mentum, and energy conservation to model the flow, 
including appropriate physical models for the high 
temperature environment. It faces most of the same 
challenges as computational fluid dynamics regarding 
cost and simplicity of generating answers to be used in 
the design process. Surface and volume grids are often 
difficult and/or time consuming to generate. Perturba- 
tions to the configuration are not easily accomodated 
unless special considerations and foresight were ap- 
plied in the initial grid generation process. Order 10 5 
or greater cells are required to discretize the flowfield 
around real configurations, and solution of the gov- 
erning equations requires hours to tens of hours per 
steady state solution on a CRAY C-90. 

Computational aerothermodynamics faces some 
challenges unique to the high temperature, hypersonic 
flow environment. Equations for chemical and ther- 
mal nonequilibrium, including source terms which may 
add stiffness to the equation set, must be included. A 
single case may include flow domains that vary from 
broad, subsonic, high pressure and temperature stag- 
nation regions behind strong bow shock waves to high 
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speed, low density (high Knudsen number) regions in 
the lee side. 

The grand challenge is to implement these solutions 
on the order of minutes so that they may be an in- 
timate and integral part of the design process. We 
seem to be a long way from this goal. Structured 
grid generation for the X-33 1 (a configuration with 
wing, tail, deflected control surfaces extending out over 
the wake, gaps, etc.) requires two to three weeks 
of effort, with close collaboration between CFD and 
grid generation specialists. Unstructured grid genera- 
tion on this same configuration requires approximately 
one to three days, but at present is only applicable 
with confidence to inviscid flow resolution. Thin-layer 
Navier-Stokes solutions for steady, laminar flow using 
a five species chemical model requires 0(10) hours to 
get decent, first cut aerodynamics on a coarsened grid. 
Converged heating levels on a fine grid generally re- 
quire a factor of three more relaxation time. We must 
do better. Some reflections are offered on the current 
design environment, CFD state-of-art, and the types 
of advances needed to achieve the target design goals. 

Historical Perspective and Technology 
Drivers 

Driving applications for external, continuum, hyper- 
sonic flow simulation in the period 1980 to 1990 in 
the United States included generic studies for aero- 
maneuvering and aeroassisted orbital transfer vehi- 
cles, 2 Space Transportation System (STS) support, 3 
defense interceptors design, 4 Pegasus, 0 and the Na- 
tional AeroSpace Plane (NASP). 6 In Europe, pro- 
grams of note during this same period include HER- 
MES, 7 Sanger, 8 and HOTOL 9 and early work on Huy- 
gen 10 probe. The Buran 11 development and flight tests 
were a major focus in the former Soviet Union. 

Driving applications for external, continuum, hyper- 
sonic flow simulation in the period 1991 to the present 
included closeouts or continuations of several of the 
programs noted above as well as X-33 1 (a technol- 
ogy demonstrator for a fully Reusable Launch Vehicle 
(RLV)), X34 (Mach 8 test vehicle), 12 Assured Grew 
Return Vehicle (ACRV), 13 Hyper X, 14 Commercial 
Experiment Transporter (COMET) 15 (later renamed 
(Multiple Experiment Transporter to Earth Orbit and 
Return (METEOR)), Mars Exploration Program, 16 
Stardust, Skipper, 17 and Pegasus XL Wing-Glove 
Experiment. 18 Recent advances in computational 
aerothermodynamics in Europe have been summarized 
by Kordulla and Morice. 19 In Japan, OREX and 
Hyflex with eventual application to the Hope Vehicle 
dominate hypersonic CFD applications. 20-23 Though 
funding is always a concern, there is certainly no short- 
age of interesting and challenging work to do. 

The period prior to 1985 was dominated by sig- 
nificant advances in space marching code capability; 
predominantly in Parabolized Navier Stokes 24 work, 


but also a significantly mature technology base in Vis- 
cous Shock-Layer methods. 25,26 This was also a period 
in which grid generation and adaption (predominantly 
structured) addressed challenges associated with flow 
simulation and resolution around complex vehicles. 
Hypersonic flows, with the presence of strong shocks 
and complex interactions presented special challenges 
that were not adequately accomodated by traditional, 
body oriented grid systems. 

The period from 1985 to 1991 was rich in algorithm 
advancement and development and/or integration of 
physical models. Major advances in upwind differ- 
encing methods and high resolution, non-oscillatory 
schemes were introduced in this period. While al- 
gorithm and physical model research continue in the 
period since 1991, these years saw much more empha- 
sis on application of existing algorithms and models 
to more complex configurations and on ways to ex- 
ploit computer architectures to accomplish these tasks. 
Algorithm innovations in this period have been sub- 
stantially driven by near term needs of programs. The 
consequences of these shifts have been both good and 
bad. An underlying theme of this paper addresses 
these consequences. 

Role of Aerothermodynamics in the 
Entry Vehicle Design Process 

The design of an entry vehicle / reusable launch ve- 
hicle such as Shuttle or planetary entry-probe involves 
the modeling and synthesis of various nonlinear cou- 
pled systems. A few important entry vehicle subsys- 
tems are propulsion, thermal protection system (TPS), 
structure/payload, avionics, and cryotanks. The inte- 
grated configuration of the necessary subsystems de- 
termines the shape of the entry vehicle. The vehicle’s 
shape and the thermal protection system are the com- 
ponents most impacted by aerothermodynamic pre- 
dictions of the entry flow environment. For a given 
shape, the aerodynamic performance characteristics of 
a vehicle (e.g. lift, drag, and static/ dynamic stability) 
are determined using various aerothermodynamic pre- 
diction techniques. The TPS, which is the interface 
between the entry flow environment and the vehicle, 
is selected, sized, and arc-jet tested based on aerother- 
modynamic predictions of quantities such as surface 
temperature, heat transfer rate, integrated heatload, 
shearloads, pressure loads, and ablation rates. 

To better understand the current role of aerother- 
modynamics in the design process, the process of de- 
signing an aeroshell for a propulsively controlled entry 
probe for landing on the surface of Mars is briefly de- 
scribed. Figure 1 is a process flow diagram for the 
design of the lander aeroshell; this diagram is repre- 
sentative of the design process used on recent entry 
probes such as Stardust or Mars/Pathfinder. The dia- 
gram shows discipline and subsystem design tasks as a 
function of time. In general, the fidelity of the analyses 
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and the fidelity of the design increases with time. A 
demarcation between Pre-Phase A and Phase A/B/C 
is noted in the figure; at this transition in the design 
process, a large increase occurs in the analysis level of 
fidelity. 

The first step in the design process is determining 
mission requirements and objectives such as landing 
a rover at a chosen location and time on the surface 
of Mars. An initial shape of the vehicle is determined 
as a function of a number of parameters such as land- 
ing accuracy requirements, payload mass and volume, 
and the launch booster payload shroud. The mass, 
aerothermodynamic performance of the shape, and the 
initial conditions as the vehicle enters the Martian 
atmosphere determine an envelope of possible entry 
trajectories. The initial aerothermodynamic database 
is generated using low-fidelity methods such as mod- 
ifed Newtonian for the aerodynamics and engineering 
correlations for the stagnation point heating. 27 Using a 
3 degree of freedom (3-DoF) trajectory analysis, an ini- 
tial design trajectory is generated to meet the mission 
requirements and satisfy subsystem constraints such 
as maximum acceleration loads or maximum temper- 
ature and heating rates for the TPS. 

In the Phase A/B/C segment of the lander design in 
Fig. 1, the various discipline and subsystem analyses 
are highly coupled and the design process is itera- 
tive. For example, propulsion requirements for control 
depend on the GN & C algorithm; the GN & C algo- 
rithm depends on knowledge of the lift and drag forces 
along the trajectory; the lift and drag forces are ob- 
tained from the aerothermodynamic predictions. The 
aerothermodynamic predictions depend on knowledge 
of the trajectory from the 6-DoF analysis. The 6-DoF 
analysis is a function of the c.g. location. The c.g. 
location depends on the packaging and mass proper- 
ties of the vehicle such as the mass of the payload, 
TPS, and propulsion system. The mass of TPS and 
propulsion system depend on the aerothermodynamic 
predictions. Thus, to determine if the various mission 
requirements are satisfied, an iterative process is re- 
quired. 

As illustrated above, aerothermodynamics is impor- 
tant in both the early and later stages of the aeroshell 
design process. In phase A/B/C, accurate predictions 
of the aerodynamics and entry heating are required 
to design the TPS and for 6-DoF trajectory analysis. 
Aerothermodynamic data fully representative of the 
flight environment are not available from ground based 
experimental facilities. Thus, high-fidelity numerical 
simulation techniques, such as Navier-Stokes for the 
continuum flow regime and Direct Simulation Monte 
Carlo (DSMC) for the rarefied flow regime, are em- 
ployed. These methods are computationally intensive 
and require accurate modeling of relevant physical pro- 
cesses to achieve good results; further, the resources 
required for a simulation usually increase as the ac- 
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Fig. 1 Design process for a Mars Precision Lander 
Aeroshell 

curacy of the prediction increases. Therefore, to limit 
costs, high-fidelity aerothermodynamic simulation is 
only employed in the later stages of the design process 
after the design has been significantly refined. 

The highly coupled and iterative nature of the entry 
vehicle design process as exemplified by the aeroshell 
for the Mars lander is a driver for aerothermodynamic 
development. The process of design involves defining 
and, then, refining and narrowing the boundaries of 
a design space until a final solution is reached. To 
achieve an optimum design, the optimum design so- 
lution must always remain within the design space 
as the design space boundaries are narrowed. Using 
current low-fidelity, engineering, aerothermodynamic 
prediction techniques, the preceding design optimiza- 
tion constraint cannot be guaranteed because the engi- 
neering methods are not accurate enough. Therefore, 
to produce optimized entry vehicle designs, new high- 
fidelity aerothermodynamic algorithms and implemen- 
tation strategies are required that are applicable to 
the earlier stages of the design process. To achieve 
this objective, the new algorithms will have to bal- 
ance accuracy, cost, and speed and take advantage of 
advanced computational platforms. 

Algorithms 

Some of the most important algorithm advances 
for the computation of hypersonic flows have been 
in the development of upwind and non-oscillatory 
schemes. 28-30 Central difference schemes with upwind 
modeled or non-oscillatory dissipation operators are 
included in this algorithm class. 31 Certainly, upwind 
schemes evolved more to address issues of tempo- 
ral accuracy and satisfy an intuitive valuation that 
consistency between physics and numerics is impor- 
tant. High resolution, non-oscillatory schemes more 
directly confront issues of accuracy in the near vicin- 
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ity of shocks and shear layers. Prior to these advances, 
flows with strong shocks ( pi/pi >~ 100) could not be 
computed accurately with shock-capturing methods; 
either the requisite dissipation excessively smeared the 
shock or Gibb’s phenomena caused oscillations to neg- 
ative temperatures on the upstream side of the shock. 
Shock-fitting methods have advanced to a lesser ex- 
tent; evolutionary environment has favored the more 
easily coded and robust shock-capturing methods. 

Much algorithm development for hypersonic flows 
has focused on treatment of source terms. A loosely 
coupled (split) formulation of the species continuity 
equations with the main equation set often allows a 
code developer to retain algorithm structure from a 
code that does not consider reacting gases. However, 
fully coupled (unsplit) approaches generally are more 
robust and converge more quickly than loosely coupled 
formulations. Within the context of fully-coupled ap- 
proaches the source term may be handled in a variety 
of ways to enhance robustness and convergence. In sit- 
uations where a pure explicit formulation is severely 
limited by time scales in one or more species source 
terms, explicit underrelaxation can be selectively ap- 
plied to species continuity equations. 32 A more robust 
approach, easily incorporated in the context of a sin- 
gle step, explicit algorithm, is to utilize point-implicit 
treatment of the source term. 33 The Lower-Upper 
Symmetric-Gauss-Seidel (LU-SGS) 34 can be modified 
to achieve good convergence; 35 though in very ener- 
getic situations problems with elemental mass conser- 
vation should be addressed. 36 Convergence rates are 
generally good on uniform meshes but have been ob- 
served to slow on highly stretched grid. All of these 
formulations sacrifice temporal accuracy for improved 
convergence with moderate overhead as compared to 
conventional, explicit formulations. 

Convergence acceleration can be addressed on both 
an algorithmic and procedural level. Examples of al- 
gorithmic approaches include various formulations of 
implicit relaxation, multigrid-methods, defect correc- 
tion, and preconditioning. Procedural approaches re- 
fer to best practice using any available algorithm and 
include practices such as mesh sequencing and solution 
sequencing. Most benefit is expected from algorithmic 
innovation 

Implicit relaxation relieves the restriction on the 
time step associated with highly stretched grids. It 
also is used to address restriction in time step asso- 
ciated with high Dahmkohler numbers in chemically 
reacting flows. Examples of implicit formulations for 
hypersonic flows in chemical nonequilibrium include 
LU, 37,38 ADI, 38 LU-SSOR, 39 as well as the LU-SGS 
methods cited previously. Convergence acceleration of 
implicit methods for hypersonic, nonequilibrium flows 
may be achieved using global preconditioning. 40 Im- 
plicit formulations often involve innovations that sacri- 
fice time accuracy for computational efficiency. In such 


cases, time accuracy may be recovered at the expense 
of sub-iterations. 41 In general, temporal accuracy in 
hypersonic, reacting flow simulations has received only 
modest attention, and non-trivial issues regarding evo- 
lution of systems with disparate time scales (convec- 
tive, dissipative, multiple reaction groups) may arise. 42 

Multigrid methods have achieved factors of approxi- 
mately 3.5 reduction in time to convergence for hyper- 
sonic flows as compared to application of an otherwise 
identical algorithm on a single, fine grid. 43 " 45 Difficul- 
ties arise from treatment of chemical source terms and 
reduction of high frequency errors in the vicinity of 
strong shocks. This performance is not as impressive 
as results obtained in lower speed domains; opportu- 
nities for significant advances probably exist. 

Local preconditioning in multidimensional flows ad- 
dress problems associated with convergence and trun- 
cation error in the limit of incompressible flow. 46 These 
problems arise, for example, in the stagnation region 
of very blunt bodies in hypersonic flow. They are 
exacerbated in the presence of significant viscous ef- 
fects. The local Mach number goes to zero in a broad 
region around the stagnation region and both trunca- 
tion error and dissipation behave poorly, particularly 
in upwind schemes. Recent advances suggest that the 
problems can be alleviated with local preconditioning, 
but much work remains to establish optimum perfor- 
mance in the vicinity of three-dimensional stagnation 
points and low Reynolds numbers. 

Computational Singular Perturbation (CSP) 47 ex- 
ploits concepts that are related to local precondition- 
ing, but focuses on optimal formulation of the chemical 
kinetic source terms. CSP is of special interest because 
the method: (1) provides physical insight regarding 
the balance of kinetic processes in hypersonic react- 
ing flows; and, (2) provides mathematical insight for 
numerical simulation of such flows to reduce compu- 
tational expense. Reaction groups are automatically 
identified which are active, in partial equilibrium, or 
frozen. Preliminary results were very encouraging; 
however, coupling the CSP algorithm to the flow solver 
with requisite modifications to the governing equation 
set was nontrivial. The concept reached acceptable 
maturity for uncoupled analyses. 48 

As noted previously, mesh sequencing and solution 
sequencing are two procedural approaches to conver- 
gence acceleration. Mesh sequencing refers to obtain- 
ing a converged solution on a sequence of finer grids, 
where each successive solution is initialized using the 
previous, coarser grid solution. Solution sequencing 
refers to initializing a simulation at one trajectory 
point using a converged solution from a neighboring 
trajectory point. Solution sequencing can be used 
in conjunction with mesh sequencing to generate a 
matrix of solutions across a trajectory for a single 
configuration. Solution sequencing should not replace 
mesh sequencing. Based on experience in running ma- 
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trices of cases for several X-33 configurations, it is 
always more efficient to restart a new solution on a 
coarse mesh initialized with a solution from another 
trajectory point (solution sequencing) and proceed to 
sequentially finer grids only after a threshhold residual 
convergence level has been achieved (mesh sequenc- 
ing). Proceeding directly to the finest grid (skip 
mesh sequencing steps) using an initialized solution 
from another trajectory point generally requires more 
computer time for a converged solution at the new 
freestream conditions. The mesh sequencing should be 
omitted in cases where only a wall boundary condition 
is changed; here, physically motivated load balanc- 
ing can be applied to concentrate CPU cycles in the 
boundary layer. 

Similar procedural strategies are applied in block, 
space marching. 49 Block marching is a natural ex- 
tension of the planar, space marching as implemented 
in Parabolized Navier-Stokes (PNS) solutions. The 
strategy is useful when (1) embedded subsonic or sep- 
arated regions may arise in domains which are other- 
wise amenable to PNS methodology; or (2) the PNS 
methodology is not sufficiently robust for the given 
problem. Downstream blocks are initialized using the 
exit plane data from the upstream blocks. This injec- 
tion speeds the development of a converged boundary 
layer profile as compared to conventional, unstaged, 
global relaxation techniques. Mesh sequencing, from 
coarse to fine, is still advised in these applications. 

Physical Models 

The choice of physical models for an aerothermodv- 
namic flow simulation is a function of the flight envi- 
ronment, TPS, and desired accuracy. The aerothermo- 
dynamic flight environment is divided by heat transfer 
mechanisms. The three primary mechanisms for trans- 
ferring energy to the surface of a hyper-velocity entry- 
vehicle are: a) energy transfer from particles (atoms, 
molecules, etc.) colliding with surface, b) energy trans- 
fer from chemical reactions on the surface (catalysis), 
and c) energy transfer via radiation from excited par- 
ticles in the flow. Mechanisms a) and b) are called 
convective heating and mechanism c) is called radia- 
tive heating. 

Accurate aerothermodynamic predictions of convec- 
tive heating require detailed knowledge of the chemical 
composition, transport mechanisms, fluid gradients, 
and TPS material response near the surface of the 
vehicle. The gas composition near the surface is a 
function of the path the particles have travelled from 
the freestream to the shock. Thus, some level of knowl- 
edge is required about the flowfield from the freestream 
through the shock to the surface of the vehicle. To 
calculate radiative heating, a detailed knowledge of 
the state of the gas throughout the flow is required. 
Examples of a few needed quantities are rotational, 
vibrational, and electronic state populations, coupling 



and excitation mechanisms between the various en- 
ergy modes, and modeling parameters for the various 
molecular and atomic radiation bands. Further, radia- 
tion transport at a given location is coupled to all the 
other locations in a line of sight. In general, radiative 
heating predictions are much more computationally in- 
tensive than predicting convective heating. 

The choice of a TPS affects the physical models 
needed for an aerothermodynamic flow simulation. 
Ablating and non-ablating materials are two impor- 
tant classes of TPS used for entry vehicle design. 
Ablating materials are generally employed on non- 
resuable planetary probes. Ablating materials min- 
imize the energy conducted into an entry probe by 
utilizing a phase change of the solid TPS material; en- 
ergy transfer from the flow to the TPS converts the 
solid TPS material to a gas which is then carried away 
by the flow. Ablation products injected into the flow- 
field and surface recession alter the flow environment. 
Thus, these processes must be modeled to obtain ac- 
curate aerothermodynamic predictions. Nonablating 
TPS is typically used on reusable vehicles such as the 
Shuttle Orbiter. Nonablating materials minimize the 
energy conducted into the vehicle by radiating energy 
away from the body and low material conductivity. 
For these materials, the flow environment is altered 
by chemical reactions that occur on the surface of the 
TPS. These surface catalysis reactions increase sur- 
face heating and alter temperature and concentration 
gradients within the boundary layer. Thus, surface 
catalysis modeling is needed to accurately predict sur- 
face heat transfer. 

Finally, the desired level of accuracy is an impor- 
tant consideration in choosing a set of physical models. 
In general, the computational cost increases as the 
accuracy of the physical modeling increases. For a 


5 of 35 


American Institute of Aeronautics and Astronautics Paper 97-2473 




given vehicle design, the level of accuracy is a func- 
tion of the available computational resources. Various 
modeling assumptions such as chemical and thermal 
equilibrium, frozen flow, tangent slab radiation trans- 
port, or a viscous shock layer are made to reduce the 
cost of a given calculation. These assumptions deter- 
mine the physical modeling requirements. 

Figure 2 is a plot of peak heating rate versus veloc- 
ity for some past, present and future missions. The 
peak nonablative heating rates and entry velocities 
range from 30 W/cm 2 at 4.5 km/s for Viking to 40,000 
W/cm 6 2 at 50 km/s for Galileo. Below 9 km/s, the pri- 
mary heating mechanism is convection. At 9 km/s 
and above, the heating is from a mixture of radia- 
tion and convection. For Galileo, 99% of the heating 
is from radiation. Most of the TPS used for these 
vehicles is ablative; SLA-561, a silicon based ablator, 
was used for Viking while carbon-phenolic was used on 
Pioneer- Venus and Galileo. A recently developed abla- 
tor, PICA, is being used for Stardust, a non-ablating 
Beryllium heatshield was used for Fire II, and non- 
ablating, silica based tiles were used on the Shuttle 
Orbiter. 

Equation of State 

In any hypersonic flow simulation, the thermody- 
namic pressure is defined as a function of species den- 
sities and temperature. The simplest formulation is for 
a perfect gas (PG), p = ( 7 — l)pe. High temperature 
effects on gas chemistry are sometimes approximated 
with an effective value of 7 which yields the correct 
density ratio across the shock . 50,51 This approach is 
useful for aerodynamic analysis if the PG option is the 
only one available in the code, but is unsatisfactory for 
evaluation of aeroheating. High temperature effects 
are more properly handled with chemical equilibrium 
(CE) or chemical nonequilibrium (CNE) models. In 
the CE models ( p = p{p,e)) the functional depen- 
dence of pressure on density and energy for a general- 
ized gas mixture (CE-G), in which the elemental mass 
fractions are known constants or are solved from el- 
emental continuity equations, can be determined by 
the method of free energy minimization or equilibrium 
constants . 52 The method of Liu and Vinokur 53 can be 
used to define thermodynamic curve fits for arbitrary 
gas mixtures based on species partition function infor- 
mation. A very helpful and recent review of partition 
models for air species has been prepared by Capitelli 
et. al 54 and Giordano et. al . 55 Curve fits for specific 
gas mixtures (e.g. air 56 (CE-air), CF 4 57 (CE-CF4)) 
when available are generally less computationally in- 
tensive. In the CNE models, pressure is defined as a 
summation of the partial pressures of each constituent, 
thermally perfect species p = J2 s p s RT/M s . The 
partial pressure of free electrons is computed using 
temperature T e if different from the heavy particle 
temperature T. The species densities are computed 


from species continuity equations. 

Chemical Kinetics 

Chemical source terms may be expressed as func- 
tions of thermal equilibrium chemical kinetic models 
(CNE-TE) or thermal nonequilibrium chemical kinetic 
models (CNE-TNE). Arhenius formulations for reac- 
tion rates are used as the basis for both CNE-TE and 
CNE-TNE models. Several models for air ° 8 and the 
Martian atmosphere 59 are available in the literature. 
Updates to the models have been reported. Some 
of these updates have been driven by results of the 
Bow Shock Ultraviolet Experiment in which measured 
radiative energy transport was in significant disagree- 
ment with current theoretical models . 60-62 In this 
regard, trace species that are important radiators may 
be omitted from the initial flowfield simulation and 
evaluated after the fact, provided their role in impor- 
tant reaction groups and the coupled radiation effects 
are negligible. 

The CNE-TNE models are most often derived from 
(CNE-TE) models by applying correction factors or 
effective temperatures to the rate constants. The cor- 
rection factors for dissociation-recombination are func- 
tions of both the vibrational and translation temper- 
atures and may include several additional parameters 
describing the molecules . 63 Correction factors for ex- 
change reactions are also available in some models. A 
very helpful and brief overview of the available models 
is contained in the recent paper by Losev 64 (21 models) 
and Kovach et. al 65 (15 models). A more comprehen- 
sive review (in Russian), associated with the AVO- 
GADRO project, is in the manual by Sergievskaya et. 
al . 66 

The CNE-TNE models need also address internal 
energy balance issues. If molecules are preferentially 
dissociated from higher vibrational energy states, the 
vibrational energy equation source term must account 
for this process. The model by Knab 67 directly ac- 
counts for this relation in the CNE-TNE model. Other 
models require empirical input. 

Thermal State of Species s 

The thermal state of species s may be modeled as 
Equilibrium (TE) or NonEquilibrium (TNE). 

The distribution of energy among translational, ro- 
tational, vibrational, and electronic modes in species 
s under equilibrium conditions is defined by the par- 
tition function as a function of a single temperature 
T. The specific heat capacity C p of species s is de- 
rived from the partition function. Curve fits of specific 
heat as a function of temperature (TE-C) for species of 
interest occurring in hypersonic flow simulations (plan- 
etary atmospheric constituents, fuels and combustion 
byproducts, ablation and pyrolisis gases) are available 
in the literature . 68 The actual partition functions may 
also be used (TE-P), particularly if they are already 
available to support thermal nonequilibrium options 


6 of 35 


American Institute of Aeronautics and Astronautics Paper 97-2473 



in the code. In continuum situations, the translational 
and rotational contributions are generally assumed to 
be fully excited and the vibrational contributions are 
modeled as anharmonic oscillators. In many instances, 
the electronic contribution to the partition function is 
omitted - only the modes up to and including vibration 
are retained (TE-Pv). This approximation is probably 
not bad for conditions when ionization and visible ra- 
diation are insignificant; however, we are not aware 
of studies defining the bounds of this approximation. 
In any event, the energy does not disappear (assum- 
ing conservation laws are formulated correctly!); it is 
simply distributed among other available modes, thus 
raising the temperature for a given energy increment. 

The TNE models may be subdivided into four sub- 
classes. 

The simplest sub-class is the two-temperature model 
(TNE2) in which it is assumed that the thermal state 
of species s and all other species in a mixture contain- 
ing some molecules may be defined by two tempera- 
tures, and effects of cross mode coupling are approxi- 
mated (e.g. rigid dumbell). The translational and ro- 
tational energy content of all heavy particles is defined 
by temperature T. The vibrational energy content of 
all molecules is defined by temperature Ty (sometimes 
called a lumped vibrational temperature model). The 
electronic and free-electron translational energy modes 
(if included at all) are usually assumed to be in equilib- 
rium with the vibrational modes at temperature Ty. 
(An exception to this approach is discussed by Olyn- 
ick, 69 in which electron translational energy modeled 
as a function of T was found to yield better agreement 
with experimental data from Project FIRE.) Justifica- 
tion for this last assumption in air is based on efficient 
energy transfer between electrons and the vibrational 
modes of molecular nitrogen. The TNE2 model may 
be implemented in two ways. If equilibrium curve fits 
are available, as discussed previously, the vibrational 
and electronic heat capacity of species s is easily de- 
fined as a function of Ty . It is obtained by subtracting 
the simple translational-rotational contribution from 
the curve fit for all energy modes appropriate for each 
atom or molecule (TNE2-C). If partition functions are 
available, the appropriate temperature is substituted 
into the function for each energy mode (TNE2-P). As 
before, if only the modes up to and including vibration 
are retained, the model is designated (TNE2-Pv). 

The next logical TNE sub-class is the three- 
temperature model (TNE3). It is identical to the 
TNE2 sub-class except that the electronic and free- 
electron translational energy partition is defined by a 
third temperature T e . (This designation is a misnomer 
if there are no molecules present in the gas.) This 
model has been considered primarily for high speed 
Earth entries characteristic of return from geosyn- 
chronous Earth orbit or above. If molecules are 
present, the use of curve fits is not possible because 


the relative contributions of vibrational and electronic 
modes to the energy content are not separable without 
some partition function information. 

The third TNE sub-class is the multi-temperature 
model (TNEM). It is identical to the TNE2 sub-class 
except that each molecular species s retains it own vi- 
brational temperature T v<s . When the thermal model 
retains only the modes up to and including vibration, 
the model is designated (TNEM-v). Such models are 
most often used in the code validation process using 
ground-based experimental data in which separate vi- 
brational temperatures have been measured. Separate 
rotational temperatures may also be considered in this 
classification. Such an approach is featured in a paper 
by Bose et. al. 60 in which the high altitude kinetics re- 
quired for accurate prediction of NO number densities 
and radiation for the Bow-Shock Ultraviolet experi- 
ment necessitated more complete characterization of 
the molecular energy content in the kinetic model. 70 
Separate electronic temperatures T e ^ s logically follow 
in this model, but we are unaware of studies involving 
(or justifying) this degree of difficulty. 

All of the TNE models to this point invoke an as- 
sumption of equilibrium partitioning of energy within 
two or more groupings of internal modes. This as- 
sumption can be relaxed further by discretizing each 
species s into multiple sub-species (excited states) as 
a function of vibrational and/or electronic quantum 
number(s). 71 In this fourth TNE sub-class, a state- 
specific approach (TNEQ), a unique translational tem- 
perature is still appropriate within the context of a 
continuum flow simulation; however, the role played 
by modal temperatures in defining energy distribution 
within the mode is replaced by a discrete account- 
ing of various excited states. Each state requires its 
own continuity equation, with appropriate models for 
transport and source terms. 

Transport Properties and Diffusion Models 

Transport properties are required to describe the 
diffusion of mass, momentum, and energy in a gas. 
A comprehensive review of the subject is provided 
in the text by Hirschfelder, Curtiss, and Bird. 72 The 
simplest definition, usually restricted to perfect gas 
simulations, employs Sutherlands law for viscosity, 
p. = ciT 3 / 2 /(c 2 +T), and Prandtl number for con- 
ductivity, k = nC p /Pr, where ci, Co, and Pr are 
constants. Thermochemical equilibrium flow simula- 
tions most often use curve fits for p, and k as a function 
of temperature and pressure. 73 Transport properties 
in thermochemical nonequilibrium flows are defined as 
functions of constituent species transport properties 
and respective mole fractions in the mixture. In- 
dividual species transport properties are defined as 
functions of collision cross sections; often this data is 
readily available in form of curve fits as functions of 
temperature. 74 Mixing laws for viscosity and thermal 
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conductivities based on work by Yos 75 and Wilke 76 are 
also discussed and compared to a simple, mole fraction 
weighted average by Palmer. 77 

The definition of diffusion coefficients should be dis- 
cussed within the context of the diffusion model. Full 
multicomponent diffusion involves expensive compu- 
tation of multicomponent diffusion coefficients as a 
function of binary diffusion coefficients. Application 
of the Stefan-Maxwell relations simplifies the requisite 
formulation of diffusion coefficients, but still involves 
significant overhead (e.g. evaluation of binary diffu- 
sion coefficients and sub-iterations) in the computation 
of multicomponent diffusion. These approaches are 
equivalent in terms of simulation accuracy and they 
rigorously conserve mass; the sum of all species conti- 
nuity equations telescopes (term by term cancellation) 
exactly to the mixture continuity equation. A sim- 
plification, originally derived and justified for trace 
species, 78 involves formulation of an effective binary 
diffusion coefficient ( I for species s using recip- 
rocal molal averages of binary diffusion coefficients 
(D St k), where the mass diffusion flux j s is driven by 
the gradient of mole fraction of species s (X s ) in the 
absence of pressure and thermal diffusion. This for- 
mulation is expressed 


js 


where 


/'4j' 1 -w s )D s , m VI s 


D, 


E 


Xk 
D s,k 


( 1 ) 

( 2 ) 


p is density, uj s is mass fraction of species s, and M 
and M s are mixture and species molecular weights, re- 
spectively. In applications involving velocities below 7 
km/s in air on the Space Shuttle, the heating com- 
ponent associated with the diffusion of atoms toward 
the surface and their catalytic recombination is ade- 
quately described by this approach. 49 However, species 
continuity equations written in this simplified manner 
do not generally telescope into the mixture continuity 
equation because the approximation to diffusive mass 
flux j s to not telescope to zero. At higher velocities 
and/or at smaller length scales where nonequilibrium 
effects associated with the catalytic recombination be- 
come more pronounced, the approximation to multi- 
component diffusion will deteriorate. 

Three approaches may be applied to correct this 
deficiency in mass conservation when approximating 
multicomponent diffusion of species s as a function 
of a single gradient in species s. The bifurcation 
approximation closely approximates actual diffusion 
coefficient data in a manner that imposes mass con- 
servation. 79 Explicit corrections on the diffusion fluxes 
may be imposed that remove a net background flux in 
Eq.‘ 1. 

js = js-^s^jk ( 3 ) 

k 





Fig. 3 Convective heating distribution on Stardust 
vehicle with various approximations to multicom- 
ponent diffusion. 


Finally, all diffusion fluxes may be reformulated as a 
function of viscosity p, and a constant Schmidt number 
Sc. 

j s = ^Xcos (4) 
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The effects of these approximations are seen in Fig. 3 
in which the convective heating to the Stardust vehicle 
are computed at an actual trajectory point (F^ = 10.9 
km/s, p = 2.6910 -4 kg/m 3 ) and at progressively lower 
velocities at the same density. The mole driven, un- 
corrected results refer to Eq. 1. The mole driven, 
flux corrected results refer to Eq. 3. At all energies, 
the flux corrected result is in excellent agreement with 
the Stefan-Maxwell relations. As flow energy decreases 
the uncorrected approximation approaches the Stefan- 
Maxwell reference, with average error decreasing from 
30% at the highest velocity, 15% at 7 km/s, and 12% 
at 6 km/s. Even at 6 km/s (Mach 19) the contri- 
bution to total heat flux from diffusion of species to 
the surface under fully catalytic conditions is quite 
large as compared to the non-catalytic surface case. 
The heating for all approximations to multicompo- 
nent diffusion are nearly identical when the diffusional 
contribution is removed at a non-catalytic wall. The 
effective Schmidt number changes across the velocity 
domain tested here, with Sc = 0.5 a reasonable ap- 
proximation at the lower velocities. 

Surface Catalysis and Ablation 

In a hypersonic flow, surface catalysis is a primary 
mechanism for transferring energy to the surface of an 
entry vehicle. Catalytic surface heat transfer occurs 
when the TPS material at the fluid-surface interface 
acts as a catalyst to an exothermic reaction involv- 
ing chemical species impinging the surface from the 
flow. A fraction of the catalytic energy released at the 
surface is conducted into the vehicle; the additional 
heating increases the TPS mass. Atomic oxygen com- 
bining to form 0-2 or atomic nitrogen combining to 
form IVo are important surface catalytic reactions in 
air. CO combining with O to form CO 2 is an impor- 
tant catalytic reaction in an atmosphere with a large 
CO 2 component such as Mars. 

Surface catalysis is an especially important heating 
mechanism for dissociated flows in chemical nonequi- 
librium. In general, the surface of a hypervelocity 
entry vehicle is much cooler than the flow tempera- 
ture directly behind the shock. For a boundary layer in 
chemical equilibrium, most of the dissociated species 
recombine before they reach the surface. Therefore, 
the majority of surface heat transfer is caused by par- 
ticle collisions with the surface (convection). For a flow 
in chemical nonequilibrium, flow species which disso- 
ciate in the shock-layer do not have time to completely 
recombine (three-body collision required for recombi- 
nation) before impinging the vehicle’s surface. Thus, 
the catalytic fraction of the total heat transfer rate is 
significant and must be considered to accurately pre- 
dict the surface heating. 

For the Shuttle Orbiter, catalytic heating is impor- 
tant because the Shuttle entry trajectory envelope is 
bounded by TPS material temperature limits on the 


nose and leading edges where the magnitude of cat- 
alytic heating is largest; 80 For RLV, this will also be 
true. Thus, the ability to predict catalytic heat trans- 
fer is directly linked to Shuttle entry performance. 
Predicted magnitudes for catalytic heating over the 
Shuttle surface are as follows. 81 Near the nose, the cat- 
alytic fraction of the heatload is 15-40% of the total 
and increases as the stagnation point is approached. 
On the wing leading edges, the catalytic fraction is 
between 10-20%. On most of the windside acreage, 
the catalytic fraction is between 5-10%. Finally, on 
most of the leeside surfaces, the catalytic fraction is 
negligible. 

For most CFD applications, the modeling of surface 
catalysis is based on simple kinetic theory. For atoms 
impacting a surface and recombining, the surface reac- 
tion rate is assumed proportional to the number flux of 
atoms striking the surface times a catalytic efficiency; 
the catalytic efficiency is the fraction of atoms that re- 
combine on the surface. These atoms are supplied by 
diffusion. Thus, in terms of the mass flux of species s, 
a catalytic surface boundary condition is given by, 82 

(5) 

where is the catalytic efficiency. A more rigorous 
approach to surface catalysis would involve modeling 
such phenomena as the availability of surface reaction 
sites, the rate at which surface reaction sites are filled 
or depleted, surface energy accommodation, and ac- 
tual surface reaction rates. 

For a number of resuable TPS materials and coat- 
ings, 7 curve-fits for atomic oxygen and nitrogen sur- 
face catalysis have been experimentally determined as 
a function of temperature 83 as, 

7* = a s exp(-|f) (6) 

Details of the experimental procedure which involves 
arc-jet testing are described by Stewart. 84 

A fully catalytic wall boundary condition is specified 
in Eq. 5 when 7* = 1 (i.e., all atoms striking the sur- 
face recombine) . A non-catalytic wall is specified when 
7s = 0 (i.e., the surface does not catalyze reactions). 
Depending on the material and the flow environment, 
the ratio of heating between a fully catalytic wall and 
a non-catalytic wall is as high as 2 or 3. 

All of the discussion on surface catalysis has as- 
sumed a reusable, nonablating TPS system. The 
physics and models for ablating and charring systems 
has recently been reviewed by Milos and Rasky 85 and 
Milos and Chen. 86 Significant, additional complexity 
can be introduced in these situations because of the 
additional chemical species introduced into the flow 
and boundary conditions dependent on a coupled ma- 
terial response. Simpler, uncoupled analyses may also 
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be accomplished in which species mass fractions and 
temperatures at a surface are supplied to the flowheld 
analysis from an external source. 

Turbulence and Transition 

Algebraic turbulence models (e.g. Baldwin- 
l.imiax. 87 Cebeci Smith* 8 ) remain the favored ap- 
proach to turbulence modeling in hypersonic flowheld 
codes. Use of local values of shear in the damping 
function formulation are generally required in these 
methods in hypersonic applications. 89 Provisions for 
separation, as in front of deflected bodyhaps, can be 
accomodated in some situations. 90 However, such 
models are not generally applicable in regions of mas- 
sive separation as occur on the leeside of vehicles at 
high angle of attack or in the wake, in corners, or in 
gaps. The algebraic models will often produce con- 
verged solutions in regions where the model is inappro- 
priate; the user must carefully interpret such results. 

One- and two-equation turbulence models 91 are 
more generally applicable. As experience and confi- 
dence increases with these models, they are replac- 
ing the algebraic models for howheld simulations over 
complex configurations. Several models used in hy- 
personic airbreathing propulsion applications in the 
LARCK code have been cited recently by Kumar et. 
al. 92 

Modeling / prediction of transition fronts has taken 
different tacks in applications associated with hyper- 
sonic airbreathing propulsion and rocket propelled ve- 
hicles. Prediction of transition is a far more critical 
issue in the case of hypersonic airbreathing propulsion 
vehicles because drag due to skin friction can severely 
limit performance, yet turbulent mixing of hydrogen 
and air in the combustor is required. State-of-the-art 
methods, recently reviewed in several sources, 93,94 are 
being developed and tested with these applications in 
mind. 95 The status of flight data to evaluate such mod- 
els was recently discussed by Bushnell. 96 

The approach to transition prediction can afford 
to be somewhat more conservative in rocket powered 
vehicles (e.g. Reusable Launch Vehicle, X-33, X-34, 
X-38) . The effect of transition on aerodynamics is less 
significant on these blunter bodies at high angle of at- 
tack. Also, peak turbulent heating levels are usually 
less than peak laminar heating levels on return from 
low-Earth-orbit. (This is not necessarily the case in 
suborbital test flights.) 

Wind tunnel data from conventional (non-quiet) 
tunnels are believed to provide conservative estimates 
for transition behavior except for cases in which ex- 
treme roughness or trips emerge on the flight hard- 
ware. Flight experience from Space Shuttle provides 
additional information on the range of conditions at 
which transition may occur on tiled surfaces and some 
instances where obstructions inadvertantly introduced 
between tile gaps caused early transition. 97 It has been 


observed that transition fronts in the wind tunnel for 
X-33 are not necessarily coincident with contours of 
momentum thickness Reynolds number. However, mo- 
mentum thickness Reynolds number appears to be a 
reasonable metric for defining the first occurrence of 
transition on the vehicle in the wind tunnel. The 
magnitude of this metric for transition prediction is a 
function of roughness height to boundary-layer thick- 
ness. Transition fronts in computation of X-33 are 
then introduced across a computational coordinate or 
zone boundary based on wind-tunnel derived metrics. 

Radiation Coupling 

A few mission scenarios where radiative heating is 
important for entry capsule design are sample returns, 
outer planet, and manned Moon or Mars missions; all 
of the mission types employ high-velocity direct entry 
trajectories for landing. Another area where radia- 
tive heating is important is rocket base heating via 
plume radiation on ascent trajectories. As an exam- 
ple, six missions are in various stages of development 
involving sample returns to Earth: Stardust, Gene- 
sis, Alladin, Mars Sample Return, Champollion, and 
Muses-C; these missions return samples from the tail 
of comet, solar wind, the moons of Mars, the surface 
of Mars, the surface of a comet, and the surface of an 
asteroid, respectively. Earth entry velocities for these 
missions vary from 12-16 km/s. For Earth entry, radi- 
ation as a heating mechanism becomes significant for 
atmospheric entry velocities above 10 km/s. 98 Thus, 
radiative heating will be important for all of these mis- 
sions. 

Computationally, radiation coupling is important 
when the fraction of flowfield energy consumed by radi- 
ation excitation is significant; the amount of radiation 
is overpredicted when radiation coupling effects are 
neglected. The radiation is overestimated because ra- 
diation coupling removes energy from the flow. If this 
energy sink is not considered in the simulation, then 
the excess energy goes toward chemical and internal 
energy mode excitation which increases the amount 
of predicted radiation. This situation is analogous to 
chemical or internal energy coupling in a flow simu- 
lation. As an example, if a perfect gas formulation is 
used to model a high-velocity entry flow, then the post 
shock temperature is grossly over-predicted. Further, 
if the perfect gas temperature is used to calculate the 
chemical composition of the gas, then the amount of 
ionization and dissociation will be overestimated. 

The amount of radiation in a flow field is a func- 
tion of the geometry, atmosphere, and entry trajectory. 
Physically, the amount of radiative heating is pro- 
portional to the shock-layer thickness which increases 
as the nose radius increases. Conversely, convective 
heating decreases as the nose radius increases; the 
magnitude of the radiative heating at the stagnation 
point is proportional to the nose radius while the con- 
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vective heating is inversely proportional to the square 
root of the nose radius." For example, Stardust (12.6 
km/s relative Earth entry velocity, 8.2° entry angle, 
60° sphere-cone, nose radius of .23 m, maximum diam- 
eter of .8 m) has a peak total non- ablative heat transfer 
rate of 1300 W/cm 2 where the radiative component is 
about 10% of the total. Fire II 100 (11.4 km/s relative 
Earth entry velocity, 14.6° entry angle, nose radius of 
.8 m, maximum diameter of .63 m), a flight exper- 
iment in early 1960’s, has a peak total non-ablative 
heat transfer of 1140 W/cm 2 where the radiative com- 
ponent is about 40% of the total. Thus, because the 
Fire II vehicle has a larger nose radius and a steeper 
entry trajectory, the magnitude of the radiative heat- 
ing and its fraction of the total is much larger for Fire 
II than for Stardust despite the lower entry velocity 
for the Fire II vehicle. 

For manned missions, vehicle dimensions are usually 
much larger than for the robotic missions. Therefore, 
convective heating is reduced, but the amount of ra- 
diative heating is increased. For example, on Apollo 
(11 km/s entry), peak heating is about 350 W/cm 2 
with half the total coming from radiation. 

Finally, the composition of the atmosphere affects 
the magnitude of radiative heating and its spectral 
distribution. For example, for an equivalent trajec- 
tory and geometry, radiative heating in air is less than 
radiative heating in a C0 2 — N 2 mixture representative 
of the Martian or Venusian atmospheres. For H 2 — H e 
mixture representative of the Jovian atmosphere, an 
entry velocity on the order of 35 km /sis needed to gen- 
erate the post shock temperatures for air at 11 km/s. 
Thus, for outer planet missions, the entry velocity at 
which radiative heating becomes important is much 
higher than for an Earth entry. 

Computationally, the high-fidelity simulation of a 
radiating flowfield is daunting; radiation modeling can 
increase the cost of a calculation by orders of magni- 
tude. At the high entry velocities where radiation is 
important, flow species are both dissociated and ion- 
ized. Thus, for a Navier-Stokes simulation with finite 
rate chemistry and thermal nonequilibrium, a large 
number of physical processes are modeled and a stiff 
equation set must be solved; this requirement is true 
regardless of whether radiation effects are included. 
With radiation, the number of physical processes that 
are modeled is orders of magnitude greater. Further, 
ablative TPS materials are typically used for vehicles 
with large entry velocities. The inclusion of ablation 
products adds more complexity because they can ab- 
sorb radiation. Thus, the radiation properties of the 
ablation products must be modeled to accurately pre- 
dict the surface radiative heating. 

Radiation intensity predictions over a relevant spec- 
tral range at a point in an entry flowfield requires 
emission and absorption coefficients for important ra- 
diation mechanisms, excited state population distri- 


butions, and radiation transport. To calculate the 
absorption and emission coefficients, a large database 
of modeling parameters is required; these parameters 
are estimated from experimental data and theoreti- 
cal formulations which are sometimes questionable for 
the velocity range or atmosphere of interest in en- 
try capsule design. To predict radiative intensities, 
excited state population distributions are needed; ra- 
diation intensity is directly proportional to the number 
of particles in an excited state. For typical gas mix- 
tures, there are thousands of excited states. If the 
excited state populations are in nonequilibrium, then 
without simplifying approximations, a rate equation 
is needed for each state. Again, this requires a vast 
database of modeling parameters. Finally, the radi- 
ation transport at thousands of spectral frequencies 
must be calculated; radiation is attenuated through 
emission and absorption as it travels from one point in 
the flow to another. Therefore, the radiation intensity 
at each point in the flow is a function of every other 
point in the flow within a line of sight; the radiation 
transport algorithm models this effect. 

In recent years, strides have been taken to im- 
prove flowfield radiation prediction and coupling tech- 
niques. For modeling nonequilibrium radiation, the 
NEQUAIR 101 and LORAN 102 codes were developed. 
NOVAR, a derivative of LORAN, has been loosely cou- 
pled with the flow solver GIANTS to generate Navier- 
Stokes solutions with coupled radiation for comparison 
with the FIRE II flight data. 69 Using a tangent slab 
radiation transport algorithm, a good agreement with 
the flight data (total stagnation heating and radiative 
intensity) over the entry trajectory was obtained; fur- 
ther, radiative coupling lead to a better comparison 
with the flight data as compared to an uncoupled simu- 
lation. Most calculations with radiation coupling have 
used a tangent-slab approximation to model the radia- 
tion transport; the assumptions of this approximation 
are valid in the stagnation region of a blunt body 
flow. Using LAURA and LORAN, however, a few 
Navier-Stokes calculations have been generated with 
coupled radiation and a multi-dimensional radiation 
transport algorithm. 103 Moving to higher levels of ap- 
proximation, flow predictions with radiation coupling 
have been performed, as non-equilibrium and 1-D 104 
and VSL with equilibrium radiation and ablation. 105 
For Hydrogen-Helium, radiation work done in the late 
1960’s would still be considered state-of-the art. 106 

The accurate modeling of radiation and radiation 
coupling is one of the great remaining challenges in 
computational aerothermodynamic design. With the 
success of Galileo, our ability to design an entry cap- 
sule and TPS to survive tremendous heatloads domi- 
nated by radiative heating was confirmed. The heat- 
shield, however, was over-designed and the pre-flight 
predictions were not very accurate based on current 
interpretation of flight data. 107 Even in cases where 
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the magnitude of radiation is low, recent flight exper- 
iments reveal serious shortcomings in some elements 
of the models. 61,62 A question that remains to be an- 
swered is: for flows with significant radiation coupling, 
can we develop prediction methodologies with enough 
accuracy to reduce design margins? To achieve this 
goal, research and development is needed to couple 
radiation and ablation to 3-D Navier-Stokes solvers, 
develop appropriate radiation process models (partic- 
ularly for other planetary atmospheres, develop multi- 
dimensional radiation transport algorithms, and con- 
tinue to improve nonequilibrium radiation modeling. 

Low Density Flows 

Computation of nozzle flows and wake flows often 
involve severe expansions into domains that are not ad- 
equately modeled with the continuum approximations 
in the Navier-Stokes Equations. Bird 108 (Eq. 8.45) de- 
fines a breakdown criteria of translational equilibrium 
in expansions based on physical arguments. Numerical 
error, particularly that associated with approximate 
Riemann solvers in these extremes, may cause the 
continuum solutions to deteriorate even earlier as man- 
ifested by entropy production. 109 

Similar problems occur in continuum simulations of 
hypersonic flows over vehicles at high altitude. Di- 
rect Simulation Monte Carlo methods offer the best 
simulation capabilities at these rarefied conditions. 
Application of models for temperature, concentration, 
velocity, and pressure slip 110 in a continuum simula- 
tion extend the altitude range for stagnation point 
heating agreement with DSMC to 115 km for the 
Shuttle Orbiter nose radius at 7.5 km/s. 111 Pitching 
moment coefficient prediction by DSMC and Navier- 
Stokes simulations for the COMET reentry vehicle at 
90 km and 60 deg angle of attack are also in very 
good agreement. 112 However, significant disagreement 
in shock thickness and shock layer profiles exist down 
to at least 90 km in these cases. 

Surface Definition and Grid Generation 

Structured Grids 

Grids are the foundation upon which all hypersonic 
CFD analysis is built and without quality grids, it is 
impossible to obtain solutions of high accuracy. The 
grid generation process represents one of the great- 
est impediments to the timely application of CFD to 
the design process. At the present, the creation of 
a high quality CFD grid is a labor intensive process 
requiring skilled workers with access to only a few au- 
tomated procedures and requiring days or weeks to 
complete. Although grid generation software has im- 
proved greatly in the past decade, it is still an “art” 
to create a quality grid. 

Surface Modeling 

All complex surfaces are generated on any num- 
ber of commercial or company proprietary CAD/CAM 


systems which theoretically should be able to commu- 
nicate through the standard IGES protocol with other 
design systems as well as grid generation software. Ex- 
perience has shown that this communication between 
CAD/CAM systems and grid generation software all 
to often is less than adequate. To alleviate this prob- 
lem, NASA created the NASA IGES 113 software that 
takes any generic IGES file and creates a file com- 
patible with grid generation software used by NASA. 
This step has improved the communication process 
but there is an additional disconnect between design- 
ers and grid generators. In general, designers using 
CAD/CAM are not sensitive to the surface quality 
requirements for grid generation. As a result, con- 
siderable time can be required to “repair” a surface 
representation before the grid generation process can 
begin. The best results are obtained when grid gen- 
erators become involved early in the design process so 
their requirements can be incorporated in the surface 
model. 

Surface Grids 

At the present time, the surface grid generation 
process is the most time consuming and labor inten- 
sive element in grid generation. It is at this point 
that the “patches” or surfaces, which can number in 
hundreds for complex configurations, must be pulled 
together, grid densities and distributions imposed on 
the grid, grid topology established, the grid smoothed 
and projected to the surface. There are a number of 
commercial, public domain and company proprietary 
software packages such as ICEM, 114 GRIDGEN, 11 ' 5 
NGP, 116 GENIE 117 and EAGLE 118 to name a few that 
can be used to generate surface grids and block faces of 
a computational domain. Most of these codes include 
a graphical interface to enhance the process. 

Parametric studies, i.e. control surface deflection, 
body shape change, etc., can require only minor 
changes in or complete rebuilds of the surface grid de- 
pending on the nature of the grid. Several of the codes 
mentioned above offer the option to build a paramet- 
ric model of the configuration and the automation of 
parametric changes in the grid at a high cost in initial 
set up time. This approach is effective when the set up 
costs can be spread out over a long term project but 
not applicable to short term “paper studies” or when 
an evolving vehicle is undergoing large and frequent 
geometric changes. 

Volume Grids 

Volume grid generation is the most automated ele- 
ment of the grid generation process. This is the point 
at which grid boundary conditions are set and the 
grid can be smoothed. This process is normally ac- 
complished in a batch mode on either a workstation 
or mainframe computer depending on the size of the 
job. The ICEM, NGP, GENIE, and EAGLE codes all 
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have a volume grid generation capability. In addition, 
GRIDGEN3D, 115 3DGRAPE, 119 3DMAGGS 120 and 
3DGRAPEAL 121 are stand alone elliptic grid gener- 
ators. Also Alter 122 has developed software that can 
be used to quickly make local alterations in existing 
volume grids. Resolution of important flow features 
which are not necessarily aligned with the surface or 
bow shock can be accommodated using solution adap- 
tive methods as represented by SAGEv2. 123 

Unstructured Grids 

While subsonic and transonic viscous analyses have 
been demonstrated with some success using unstruc- 
tured grids, such efforts for hypersonic flows with heat 
transfer are less evident. 124 There is, however, a ma- 
turing technology base for application of unstructured 
grids to inviscid, hypersonic flow analyses. 

Although unstructured grid generation suffers from 
the same surface definition problems encountered in 
structured grid generation, once the initial set up is 
done the construction of subsequent surface and vol- 
ume grids is quick and efficient - on the order of hours 
- which makes this approach attractive for parametric 
analysis. There are two basic types of unstructured 
grids, the traditional based on tetrahedra and carte- 
sian based on hexahedra. 

Surface Grids 

Like the structured grid process, connectivity must 
be established between surface patches and grid distri- 
butions specified before the surface grid can be gener- 
ated. Most software has an interactive front end to aid 
in the process. In general, this element of the process is 
much less labor intensive and more automated relative 
to structured grid generation. The FELISA 120 and 
VGRID 126 codes represent mature software for tradi- 
tional unstructured grid generation and the algorithm 
described by Aftosmis et. al. 127 is representative of 
unstructured cartesian software. 

Volume Grids 

Once the boundaries of the computational domain 
are established and parameters are set to control grid 
spacing, the grid is created in a batch mode. FELISA, 
VGRID and unstructured Cartesian solvers all possess 
self-consistent volume grid generators. 

Hardware Requiremments 

Except for the most complex grids, a high end work- 
station with 750 Mbytes of memory is sufficient for 
structured and unstructured grid generation. For ex- 
ample, a 127 block volume grid having a total of 5 
million grid points was generated on a Silicon Graph- 
ics /R10000 machine with 750 Mbytes of memory. 


Validation and Error Estimation 

General Comments 

Validation and error estimation are critical chal- 
lenges within computational aerothermodynamics be- 
cause uncertainties in predicting vehicle aerothermo- 
dynamic performance increase design margins; large 
design margins add weight and lower vehicle perfor- 
mance. 128 Validation requires the user check that (1) 
the governing equations are being solved correctly and 
(2) the correct and sufficient set of equations (physical 
models) are being solved. The first set of checks in- 
volve things like: temporal and/or spacial grid conver- 
gence studies; conservation checks on elemental mass 
fractions, entropy, and total enthalpy; comparisons to 
known analytic solutions for simple, well-posed prob- 
lems; independent recoding of program modules; and 
comparisons to experimental data with insignificant 
uncertainty in physical models. This aspect of vali- 
dation has much in common with other flow domains, 
has been described in the literature, 129 and will not be 
elaborated on here. The second set of checks predom- 
inantly involve comparisons to experimental data in 
which effects of empirically or approximately modeled 
physical phenomena are manifest in the measurements. 
The requisite physical properties and models for com- 
putational aerothermodynamics have been discussed 
in a previous section. Uncertainties in these physical 
properties and models tend to escalate with increasing 
energy of the flow. The difficulties and cost of obtain- 
ing experimental data for code validation also escalate 
with increasing energy of the flow. 

Error estimation is predominantly based on code 
validation experience. In any given design applica- 
tion, the code is validated against experimental data 
sets that most closely match configuration and flight 
parameter space (Mach number, Reynolds number, to- 
tal enthalpy and pressure, wall temperature, surface 
roughness, etc.). Under the most ideal circumstances, 
error estimation is based on fully grid converged so- 
lutions in which the application parameter space is 
fully contained within a validation test set parameter 
space. These ideal circumstances are often attainable 
in CFD applications ranging from incompressible to 
supersonic. At hypersonic conditions, extrapolation 
of error estimate beyond any available validation set 
parameter space is usually required. 

Grid convergence studies are a necessary but insuf- 
ficient metric for establishing error estimates in the 
less than ideal circumstances that usually prevail in 
hypersonic applications. In geometrically simple but 
physically challenging validation experiments grid con- 
verged solutions may be realized and attention can be 
focused on the physics. Experimental programs 130, 131 
and workshops have been initiated to address these 
component focused needs while trying to minimize the 
possibility that coupled modeling errors may fortu- 
itously cancel. Even with simple shapes, such studies 
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are not trivial when the entry velocities (energies) are 
sufficiently high (e.g. Galileo, Stardust, Apollo) that 
radiation and ablation must be included in the simu- 
lation. 

In more geometrically complex configurations, ob- 
taining a grid converged solution is either not possible 
or is precluded by higher priority claims on computer 
resources. (This situation is changing with advances 
in workstation clusters and validation data sets like 
the STS OEX data will need to be revisited.) In 
these instances the validation data set used to establish 
the uncertainty in computed heating levels or aerody- 
namics should be defined along with an appropriate 
description of the relevant parameter space and level 
of grid convergence. The sensitivity of the result to 
model parameters can be defined. Evaluation of the 
relative contributions to heating rate (thermal con- 
duction, diffusion, radiation) or aerodynamic forces 
(pressure, skin friction, base pressure) can sometimes 
provide insight into the error estimation process. 

Ultimately, confidence in CFD predictions depends 
on comprehensive comparisons to experimental data 
with well defined error bounds. The greater the in- 
ventory of validation sets, the better one is able to 
deal with complexity and to make error estimates. 
Consider, for example, the large flight aerothermal 
database that exists for the Shuttle Orbiter. At what 
point in a series of validation tests using this data is 
a code validated for aerodynamics and heating predic- 
tions? Certainly a single test case, even one with a 
single grid refinement that produces changes in aero- 
dynamics and heating less than experimental uncer- 
tainty, is only a starting point. Isolating physical 
modeling errors from discretization errors in a single 
test case is difficult because of the complexity, size, and 
large number of coupled physical processes needed to 
simulate the flowfield around the Shuttle Orbiter. Dis- 
cretization and physical modeling errors are coupled 
because the physical modeling is a function of flow pa- 
rameters (velocity, temperature, pressure, etc.) and 
their gradients which are a function of grid resolution. 
Further, many of the physical processes considered 
such as vibrational excitation and dissociation are cou- 
pled. Thus, errors from one physical model propogate 
to other models and isolating errors is difficult; good 
agreement with experiment may simply represent for- 
tuitous cancelling of errors. 

Taking the question a step further: Are consistent 
comparisons with experimental data along multiple 
points along a single trajectory sufficient for valida- 
tion? Are multiple points along several different flight 
trajectories sufficient? The multiple tests now begin 
to provide a statistically meaningful basis for error 
estimation for a particular configuration and parame- 
ter space. Fortuitously cancelling errors are less likely 
though still possible. 

What happens if a new configuration is proposed in 


which bodvflap deflections are larger than the Shuttle 
Orbiter at hypersonic, reacting gas conditions? If one 
has high enthalpy wedge experiments in the validation 
inventory, then reasonable error estimates may be de- 
rived. (Note that data sets for wedge experiments may 
need to account for three-dimensional effects. Also, 
some recent experimental data has suggested more 
uncertainty in established models for what had been 
thought to be a “simple” problem. 132 ) 

Code validation is a process. The process continues 
until no statistically significant reduction of error es- 
timate can be realized for a specific application. Even 
in the limit of infinitely fine grid, the error estimate 
cannot go to zero if experimental measurements are 
required to anchor physical models within the simula- 
tion. Simulation uncertainties depend on the quantity 
being predicted (aerodynamics, heating rate, electron 
number density) and increase with the complexity of 
the flow. The complexity of the flow is a function 
of both geometry and entry velocity. For example, 
simulations involving steady, laminar flow of a per- 
fect gas or a gas in thermochemical equilibrium over 
a blunt, spherically capped cone have the greatest 
confidence levels, with uncertainties predominantly as- 
sociated with discretization errors. For such cases in 
the supersonic to hypersonic regime, aerodynamic co- 
efficients would generally be expected to be well within 
two to three percent. The uncertainty errors associ- 
ated with heating rate predictions in a high energy 
flow in thermal and chemical nonequilibrium for the 
same geometry will be much higher; but no blanket 
statement can be offered for all codes under all condi- 
tions. 

Margins - A Sample Case 

Figure 4 is a plot of heating margin estimates for 
the Stardust sample return capsule. These margin 
estimates are based on fully catalytic, non-ablating, 
axisymmetric Navier-Stokes calculations performed for 
the design of the capsule. For the forebody, the major 
sources of margin are from the computation (physi- 
cal modeling and discretization errors) and angle-of- 
attack effects. A 25% computational margin is esti- 
mated for the forebody and 10% angle-of-attaek effect 
is estimated for the shoulder. These margins are cu- 
mulative. Thus, at the shoulder, the adjusted heat- 
ing value is 1.38 times larger than the zero margin 
value. In the afterbody region, the uncertainties are 
much higher. Unknowns related to turbulence, abla- 
tion, angle-of-attack, and the laminar wake heating 
computations are all major contributors to the design 
margins. At the PICA/SLA seal which is the joint 
between the forebodv and afterbody TPS, the design 
margin is a factor of 4.5. 

To estimate and develop an acceptable risk assess- 
ment for an entry vehicle requires design margins. For 
example, a TPS designer wants to know the error 


14 of 35 


American Institute of Aeronautics and Astronautics Paper 97-2473 




Fig. 4 Approximate entry heating design margins 
for the Stardust Sample Return Capsule 

bars on an entry heating calculation to determine TPS 
mass margins. At present, a formal methodology does 
not exist for determining computational aerothermo- 
dynamic design margins. The current state-of-the-art 
consists of general estimates based on comparisons 
with flight data, experimental data, previous design 
experience, atmospheric uncertainties, and limited pa- 
rameter studies. This “ad hoc” approach leading to 
the estimates in Fig. 4 is very limited. 

A formal heating margin analyis methodology is 
possible using the approach employed in trajectory 
analysis. In trajectory analysis, a state vector, con- 
sisting of important simulation parameters, is selected. 
Next, uncertainties in the state parameters about a 
marginal value are estimated. Finally, perturbed com- 
binations of the state vector parameters are selected 
via a Monte Carlo approach and a large number of 
trajectory analyses are performed. Variations about 
the reference trajectory determine design margins. For 
a Navier-Stokes heating calculation, parameters relat- 
ing to the physical models such as chemical rates, 
transport coefficients, Schmidt number, and surface 
catalycity could be selected and their associated un- 
certainties could be estimated. Then, combinations of 
these parameters would be selected followed by a large 
number of Navier-Stokes calculations. For a simple ge- 
ometry such as the forebody of a sphere and a limited 
chemistry set such as five species air, this type of anal- 
ysis, while computationally intensive, is possible with 
present day computers. 

Codes 

The following codes are representative of capabil- 
ities available for hypersonic flow analyses, includ- 
ing options for thermochemical equilibrium and non- 
equilibrium. It is not an all inclusive list; rather, it 
is restricted to codes with which we have application 
experience. 


LAURA (The Langley Aerothermodynamic 
Upwind Relaxation Algorithm) 

Inviscid flux definition in LAURA employs Roe’s 
averaging 29 and Yee’s Symmetric Total Variation Di- 
minishing (STVD) 133 for second-order accuracy away 
from discontinuities. Harten’s entropy fix (eigenvalue 
limiter) is applied across cell faces. Special variations 
of the limiter are employed across viscous dominated 
boundary and shear layers. These treatments over- 
come problems often encountered with the baseline 
Roe’s method regarding the “carbuncle” phenomenon 
or the baseline Harten’s method in which numerical 
dissipation (proportional to a large eigenvalue) com- 
petes with physical dissipation. Central differences 
are used to define viscous flux. Point-implicit relax- 
ation of the steady form of the conservation laws at 
each control volume in a computational plane is im- 
plemented; sweeping from plane to plane in a block 
and from block to block across the entire domain of 
interest. Consequently, the relaxation algorithm is Ja- 
cobi within a plane, and Gauss-Seidel from plane to 
plane. The STVD formulation of the anti-dissipative, 
second-order corrections are uniquely well suited to 
point-implicit relaxation of the steady equations; other 
formulations have been found to be less robust in 
this context. Evaluation of the inviscid Jacobian with 
eigenvalue limiter parameter e set to 0.5 increases ro- 
bustness of the relaxation scheme without detriment 
to solution accuracy. 

The LAURA code is unique among hypersonic flow 
solvers in two regards. First, it is the only known 
flow solver to exploit macrotasking on Cray computers. 
Macrotasking refers to parallel execution of code on a 
shared memory machine on a subroutine level. (In con- 
trast, microtasking and autotasking generally occur on 
a do loop level within a single subroutine.) The use of 
macrotasking, along with provisions for asynchronous 
relaxation, enable exceptionally high average concur- 
rency of tasks for a multitasked job. Asynchronous 
relaxation also enables physically motivated load bal- 
ancing in which CPU cycles are concentrated in regions 
known to converge more slowly (e.g. separated flow 
regions preceding deflected body flaps, near wake). 
Second, a built in grid adaption routine simultaneously 
aligns the outer boundary of the computational grid 
with the captured bow shock and enforces near wall 
grid resolution required for aeroheating analyses. This 
feature greatly simplifies the grid generation process 
for multiple cases over a wide range of Mach numbers, 
Reynolds Numbers, and angles of attack on a single 
configuration. 

The basic features of LAURA. 4.1 (214 page user 
manual 134 ) include options for five thermochemical ki- 
netic models for 11 species air, two equilibrium air 
models , 53,56 two-temperature thermal model (TNE2- 
C), two algebraic turbulence models, six models for 
wall catalysis, a radiative equilibrium wall, discretiza- 
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tion on up to six, simply connected, structured blocks, 
and mesh sequencing. Solutions of the Navier-Stokes, 
Thin-Layer Navier-Stokes, and Euler equations can 
be implemented. Solution Jacobians may be stored 
out-of-core for significant reduction (as much as factor 
2 - function of number of unknowns) in memory re- 
quirement. New features offered in LAURA. 4. 4 (user 
manual not yet updated) include options for Martian 
atmospheric chemistry, 1000 structured blocks with in- 
teger stride connectivity (175 block solution tested), 
and post-processing files generated to support integral- 
boundary-layer analyses of variations in surface catal- 
ysis and emissivitv of the thermal protection system. 
A PVM/MPI version of the code is being tested which 
currently supports all of the options in LAURA. 4.1 
and most of the options in LAURA. 4. 4. 

Multigrid relaxation had been tested on an earlier 
version of LAURA. 44 In some cases, it yielded con- 
verged solutions in 1/3 of the computational time at 
the expense of increased storage. However, this option 
is not currently supported because of memory over- 
head requirements. 

Code validation runs for LAURA have been docu- 
mented with STS flight data 49, 135 as well as with other 
flight test and ground based data. 136,137 

GASP (The General Aerodynamic Simulation 
Program) 

Inviscid flux definition in GASP employs several op- 
tions, including Roe’s and Van Leer’s upwind biased 
formulations and central differencing with artificial 
viscosity. Central differences are used to define vis- 
cous flux. A rich variety of options are supported 
for convergence of both steady and unsteady Navier- 
Stokes, Thin-Laver Navier-Stokes, Parabolized Navier- 
Stokes, and Euler equations. These options include 
mesh sequencing, preconditioning, approximate fac- 
torization (AF), Line Gauss Seidel (LGS), General- 
ized Minimal Residual (GMRES), 138 mesh sequenc- 
ing and multi-grid. Both algebraic and two-equation 
turbulence models with wall function options are sup- 
ported. Generalized Zonal-Boundary Interpolation is 
supported across zonal intersections defined by a sin- 
gle logical boundary. Parallel processing is supported 
on shared memory computer architectures. A compre- 
hensive set of thermochemical kinetic models is offered 
for air chemistry, hydrogen-air combustion, and hydro- 
carbons in a database containing 455 reactions and 
34 species. Thermal nonequilibrium may be mod- 
eled using a separate vibrational temperature for each 
molecule (TNEM-V) or a lumped vibrational tem- 
perature common to all molecules (TNE2-Pv). The 
comprehensive GASP V3 Users Manual 139 (654 pages) 
and the Graphical User Interface (GUI) for problem 
setup and data manipulation make GASP more user 
friendly. 

The GASP code has been validated at hyper- 


sonic, nonequilibrium flow conditions with STS heat- 
ing data. 140 

GIANTS (Gauss-Seidel Implicit 
Aerothermodynamic Navier-Stokes code with 
Thermochemical Surface conditions) 

The GIANTS code 69 is a laminar, 2- 
D/axisymmetric, full Navier-Stokes code which 
simulates thermal and chemical nonequilibrium. It 
was developed to accurately simulate the flowfield 
and surface boundary conditions characteristic of a 
pyrolizing or ablating TPS. Loosely coupled with a 
material response code, it has been used to perform 
aerothermodynamic analysis and TPS sizing for Mars 
and Earth entries; GIANTS was the primary TPS 
analysis code for the Mars/Pathfinder and Stardust 
heatshield designs. GIANTS has been coupled 
with a radiation module, NOVAR (Nonequilibrium 
Optimized VectorizeAble Radiation process model), 
and validated against data from the Fire II flight 
experiment. Also, code to code comparisons between 
GASP, LAURA, and GIANTS have been performed; 
for similar flow modeling (transport, chemistry, etc.), 
the codes agree well. 

GIANTS employs a Gauss-Seidel line relaxation 
technique based on the work of Candler and MacCor- 
mack 141 and is a derivative of a code developed by 
Candler. It employs species density equations and fi- 
nite rate chemistry to model chemical nonequilibrium, 
a vibrational energy equation to simulate thermal 
nonequilibrium, a bifurcation formulation for species 
diffusion, a modified Steger- Warming inviscid flux for- 
mulation and central differences for the viscous fluxes. 
The method is extremely efficient and robust for cal- 
culating hypersonic blunt-body flowflields. Efficiency 
and robustness are important because the integrated 
heatload is needed to size the TPS; therefore, heat 
transfer predictions are required at many points along 
an entry vehicle’s trajectory. 

DPLUR (Data-Parallel Lower-Upper Relaxation) 

DPLUR 142 is based on a variation of the LU-SGS 34 
algorithm that is designed to optimize performance on 
a single-instruction, multiple-data (SIMD) computer 
architecture. The program seeks to exploit built- 
in, optimized interprocessor communications paths by 
appropriate assignment of data to processors and re- 
quiring only nearest neighbor communication. The 
method has been applied to inviscid, single block rep- 
resentations of the forebodv of the X33 as well as 
other geometrically simpler but physically more com- 
plex simulations. 60,132 Some aspects of the method 
may not translate well into unstructured grids because 
of the architectural constraints on nearest neighbor 
mapping. 
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FELISA_HYP 

The FELISA_HYP 143 code employs an unstructured 
grid algorithm specifically constructed for robust, hy- 
personic flow simulation. It is a finite-volume based 
formulation that employs an efficient edge data struc- 
ture. Second-order accuracy is maintained in smooth 
regions using linear reconstruction following MUSCL 
concepts. 28, 144 The Local Extremum Diminishing cri- 
teria 145 is used as a limiter near flow discontinuities. A 
simple, forward Euler explicit time stepping is used to 
relax the equations. The code is currently limited to 
inviscid flows. Options for equilibrium air chemistry 
are available. 

In spite of its current limitation to inviscid flows, 
FELISA_HYP has proven particularly valuable in the 
CFD design environment because of the relatively 
quick grid generation capability. In the X-33 design 
Phase I, the FELISA grid tools could be applied by an 
experienced user to generate unstructured surface and 
volume grids in four days (first configuration took 1.5 
weeks). In comparison, a multiple-block structured 
grid for LAURA took approximately four weeks to 
construct using state-of-the-art grid generation tools 
and expertise (first configuration took 6 weeks). (The 
LAURA grid required boundary-layer resolution for 
viscous flow, but this additional constraint is thought 
to add less than 30% penalty to the grid generation 
costs.) The actual FELISA_HYP solver is somewhat 
slower than the inviscid version of LAURA; however, 
the FELISA_HYP solver would finish several inviscid 
solutions on a new configuration before LAURA could 
even get started with a usable grid. 

OVERFLOW 

The OVERFLOW code is not specifically con- 
structed for hypersonic flow simulations. The code 
is unique, however, in its use of overset structured 
grids to discretize complex configurations. This de- 
gree of freedom should decrease the grid generation 
time for entry vehicle design, provided that no special 
provisions are required for resolution of strong shocks 
traversing the overset grid regions. OVERFLOW has 
been used in high supersonic freestream conditions. 146 
The chemistry model does not now include equilibrium 
or nonequilibrium option, but does allow for frozen 
flow of multiple species (single phase), loosely coupled 
to the flow equations, to approximately account for 
solid rocket booster (SRB) plume effects on Shuttle 
ascent. The flow equations have a variable gamma 
(mixture of perfect gases/frozen chemistry) capability. 
The hypersonic simulation capability is not very ma- 
ture at present. 

3D Non-Body-Fitted Cartesian Euler 

The 3D Non-Body-Fitted Cartesian Euler flow 
solver 147 (see Gridding Strategies), like OVERFLOW, 
is not specifically constructed for hypersonic flow simu- 


lations. The flow relaxation algorithm uses a Jameson- 
Schmitt-Turkel scheme 148 using multi-stage Runge- 
Kutta integration and central differencing with added 
second- and fourth-order artificial viscosity. The 
scheme employs a feature detection algorithm which 
provides for cell enrichment across strong shocks. Al- 
though the code is currently restricted to perfect-gas 
simulation, it has been applied to a Mach 10 flow sim- 
ulation over an X-33 configuration. Computed results 
for aerodynamics appeared to be of acceptable quality 
for preliminary design study; extensive validation and 
grid refinement tests are still required. It appears to 
be competitive with the FELISA_HYP unstructured 
approach; however, more comprehensive studies are 
required. 

We believe it is worthwhile to bring codes like 
OVERFLOW and the Cartesian Euler to bear on 
hypersonic simulation problems so that the tradeoffs 
on various approaches to flow discretization (unstruc- 
tured, patched structured, overset structured, Carte- 
sian unstructured) can be understood and perhaps rec- 
ognize synergistic, hybrid approaches that best meet 
design needs. 

Applications 

For simple geometries (e.g. capsules, planetary 
probes) CFD methods can generally define surface and 
volume grids and generate a matrix of solutions faster 
than experimental, ground-based methods. Support 
for COMET 112 confirms this observation. As geome- 
tries grow in complexity, the situation reverses. At 
present, ground-based experimental methods can con- 
struct and test relatively complex configurations across 
a matrix of Mach numbers and Reynolds numbers 
much faster than structured, viscous CFD simulation 
methods. Unstructured, inviscid simulation matrices 
can be filled on a time frame that is only begin- 
ning to be competitive with ground-based methods. 
(Ultimately, if a fine-grained test matrix is required, 
involving sweeps in pitch, yaw, roll, and control sur- 
face deflections for a set configuration, even the un- 
structured, inviscid methods are not competitive with 
ground-based methods. Of course, the converse ob- 
servations can be made for simulation matrices that 
extend beyond the Mach - Reynolds number capabil- 
ities of ground-based facilities.) Experience in Phase 
I and II for X-33 confirms these observations. Nei- 
ther ground-based nor computational analyses have a 
lock on fidelity of simulation in the hypersonic flight 
environment. 

The applications that follow describe the role com- 
putational aerothermodynamics has played in defin- 
ing the aerothermodynamic environments for a variety 
of configurations and physical modeling requirements 
ranging from simple to complex. The examples pre- 
dominantly reflect the authors’ personal applications 
experience. 
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COMET/METEOR 

Aerodynamics and surface heating for the Com- 
mercial Experiment Transporter (COMET) at several 
points along its trajectory on return from Low Earth 
Orbit were calculated with LAURA and a Direct Simu- 
lation Monte Carlo method (DSMC). 112 The COMET 
module (later renamed METEOR) has no active con- 
trol system, and relies entirely on aerodynamic forces 
for stability and proper orientation during its maxi- 
mum heating pulse. The aerodynamic data base was 
used within a six degree-of-freedom trajectory code to 
define a splashdown footprint. A Direct Simulation 
Monte Carlo method was used to define the flowfield 
in a transitional, rarefied regime (above 90 km.) The 
LAURA program was used to define the flowfield in 
the transitional to continuum regime (below 90 km.) 
Wake flows were included prior to the peak heating 
point because of the large initial angle of attack. Con- 
tinuum and rarefied aerodynamic predictions for lift, 
drag, and moment were in good agreement at 90 km. 
Thermochemical nonequilibrium models including 7 
species for air were used down to Mach 15. Both vis- 
cous and inviscid solutions were used below Mach 15. 
Wake flows were included at Mach 1.5 to account for 
important base flow effects on aerodynamics. A ma- 
trix of 46 solutions was completed between February 
14 and March 23, 1995. This matrix included 10 re- 
acting, viscous flows with wake; 13 reacting, viscous 
flows without wake, 6 perfect gas, viscous flows with 
wake, and 17 perfect gas, inviscid flows without wake. 
Angles of attack varied from 0 to 90 degrees. A solu- 
tion adaptive grid was employed to swing the extended 
grid in the wake around the body behind the base at 
zero degrees to off the side at 90 degrees. Maximum 
job size was 71.2 MW on the C-90 and required 8.8 
hours on a 72 x 36 x 64 cell domain. Actual time 
on the computer for this case was only 0.98 hours be- 
cause of extensive use of asynchronous macrotasking 
relaxation. The large average concurrent CPU usage 
enabled fast turnaround for this large matrix of cases. 

Computed results were obtained prior to initiation 
of the wind tunnel test program and were in excellent 
agreement with wind tunnel data at Mach 6 (Fig. 5). 
Mach number independence for blunt body aerody- 
namics with minimal high temperature gas effects is 
evident in these results. Flight data is not available 
because the mission had to be aborted on ascent. 

Mars Pathfinder 

CFD analysis of the Mars Pathfinder was the pre- 
dominant source for aerodynamic coefficients and 
heating environment. 27,51 The spherically blunted, 70 
degree half angle cone shape is very similar to the ear- 
lier Viking probe. While a significant data base from 
the Viking Project was used, the entry parameters for 
Pathfinder (velocity and trim angle-of-attack) 149 differ 
significantly from Viking. 
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Fig. 5 Computed aerodynamic coefficients for the 
COMET re-entry vehicle at the Mach 5 and 10 
trajectory points are compared to ground-based 
experimental data at Mach 6. 

This simple shape showed significant gas chemistry 
effects on aerodynamics. As the Mars Pathfinder 
probe descends through the Martian atmosphere the 
minimum value of the post-shock effective 7 first de- 
creases from frozen gas chemistry values (« 1.333 - 
assumed effective 6 degrees of freedom in linear tri- 
atom) to equilibrium values (1.094) corresponding to a 
velocity of 4.86 km/s. As the probe continues to decel- 
erate through an equilibrium post-shock gas chemistry 
regime, the value of 7 increases again, until reaching 
its perfect gas value of 1.333 at parachute deployment 
(0.42 km/s). At small angles of attack (a < 5°) the 
sonic line location shifts from the shoulder to the nose 
cap and back again on the leeside symmetry plane be- 
cause of the change in 7 and the cone half-angle of 
70°. At 2° angle of attack, the flat, leeside pressures 
approaching the shoulder (when the sonic line sits over 
the nose) can exceed the rounded windside pressures 
approaching the shoulder (when the sonic line sits over 
the shoulder). The net effect of this crossover dis- 
tribution near the shoulder tends to pitch the probe 
to higher angles of attack. Conditions for a positive, 
destabilizing moment coefficient derivative C m _ a occur 
twice in the Mars Pathfinder mission (see Fig. 6) as 
determined by the viscous, thermochemical nonequi- 
librium flow simulations used in this study. The first 
occurrence (7.5 > LA > 6.5 km/s, 51 > h > 37 
km, vicinity of peak heating for this trajectory) results 
from the transition in the sonic line location as a func- 
tion of gas chemistry changing from nonequilibrium to 
equilibrium. The second occurrence (4.0 > LA > 3.1 
km/s, 25 > h > 22 km) results from the transition in 
the sonic line location as a function of decreasing flow 
enthalpy in an equilibrium gas chemistry regime. 

Sonic line movement also affects the heating distri- 
butions by altering the effective radius of curvature 
of the body. Peak heating tends to decrease as an- 
gle of attack increases for fixed freestream conditions 
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Fig. 6 Pitching moment coefficient versus veloc- 
ity at several angles of attack a for the Pathfinder 
entry into the Martian atmosphere. 

on the Mars Pathfinder probe because the body ap- 
pears blunter to the oncoming flow. However, the drag 
coefficient decreases with increasing angle of attack, 
so that the ballistic coefficient increases and the peak 
heating point along the trajectory may be more severe. 

Stardust 

Stardust is a Discovery class mission that is sched- 
uled for launch in early 1999 and return to Earth 
in early 2006. The prime scientific objective of the 
mission is to rendezvous with Comet Wild- 2 , gather 
samples of cometary particles and return them to 
Earth for further analysis. The captured cometary 
particles are stored in the Stardust sample return cap- 
sule. The sample return capsule enters the Earth’s 
atmosphere at approximately 13.5 kilometers/sec and 
decelerates to .6 kilometers/sec in two minutes. At 3 
km, a parachute is deployed for a land based recovery 
in Utah. Significant milestones for the mission are the 
first attempt at a sample return beyond the Moon and 
the fastest Earth entry of a human-made object ever 
attempted. Further, a revolutionary new light-weight 
heatshield material, PICA, is being flight qualified on 
the Stardust mission. 

The Stardust mission poses a number of aerother- 
modynamic challenges. Accurate predictions of the 
entry aerothermodynamic environment are needed for 
the following tasks: 1) forebody 100 and afterbody TPS 
sizing, 151 2) arc-jet testing for flight qualification of 
the heatshield, 152 ’ 153 3) predicting the aerodynamic 
stability of the vehicle during entry, 154, 155 and 4) es- 
timating the landing footprint for vehicle recovery. 100 
The flow environment for the high speed Stardust en- 
try is complex; the sample return capsule heatshield 
is ablating and the flow is radiating. Further, shape 
change reduces the drag and must be considered for 


Stardust Sample Return Capsule Design 

• Launch 1999 

• Rendevous with Comet Wild-2 



PICA/SLA Seal 


Fig. 7 Temperature contours at peak heating for 
the Stardust sample return capsule 

the footprint calculations. For the Stardust sample 
return capsule design, these aerothermodynamic chal- 
lenges are addressed using high-fidelity CFD analysis. 

Aerothermo dynamics estimates for the first two 
tasks, TPS sizing and arc-jet testing, are generated 
using the GIANTS flow solver, the FIAT 157 mate- 
rial response code, and the NOVAR radiation code. 
These modules are loosely coupled to perform flow- 
field calculations with radiation and ablation for both 
the forebody and afterbody heatshields. Fig. 7 is a 
plot of temperature contours around the vehicle for 
a typical flow calculation. For the forebody PICA 
TPS sizing, an ablation and recession model is de- 
veloped from arc-jet test results. The FIAT ablation 
model parameters are selected to reproduce recession 
and thermocouple data from a PICA arc-jet test series. 
For a given heat transfer rate, FIAT predicts surface 
temperature, surface mass fractions, and mass blow- 
ing rate; these values are used as GIANTS inputs. For 
the GIANTS flow simulation, ablation products from 
the forebody are allowed to propogate into the wake 
region; the afterbody contains a vent and contamina- 
tion from the forebody ablation products is a concern. 
The afterbody heatshield is composed of SLA -561; 
mass injection from this material is not considered in 
the simulation. 

Full-body flow calcuations are generated along the 
trajectory to define the heat pulse and size the TPS; 
the predicted thickness and total surface recession of 
the PICA heatshield are about 5 cm and 1 cm at the 
stagnation point. The effects of shape change are not 
considered in the flow simulation because the total re- 
cession at the stagnation point is only 4% of the nose 
radius. The drag reduction, however, from the shape 
change is considered in the footprint calculation. 

For arc-jet testing and heatshield qualification, es- 
timates of the non-ablative heat transfer are required; 
predicted values at a few vehicle locations are shown 
in Figs. 4 and 7. The peak non-ablative heat transfer 
at the nose assuming a fully catalytic surface is about 


19 of 35 


American Institute of Aeronautics and Astronautics Paper 97-2473 



1300 W/cm 2 with radiation. CFD predicted surface 
heat transfer, shear, and pressure values are used to set 
arc-jet test conditions. Using these conditions, arc-jet 
models are tested to estimate the performance of var- 
ious components of the heatshield design such as the 
seal between the forebody and afterbody heatshields, 
the parachute seal, the afterbody vent, and the fore- 
body TPS at stagnation conditions. 

Aerodynamic stability and the footprint calcula- 
tion (items 3 and 4) are critical to the success of 
the Stardust mission. These issues are addressed us- 
ing a 6-Degree of freedom trajectory analysis code 
(POST). The accuracy of the POST simulation and 
dispersion analysis is a direct function of the fi- 
delity of the aerodynamic database from initial at- 
mospheric entry through parachute deployment. The 
flow regimes spanned by the aerodynamic database 
are free-molecular, transition, hypersonic, supersonic, 
transonic, and subsonic. Thus, a wide variety of 
simulation tools and resources are employed. Direct 
Simulation Monte Carlo (DSMC) is used in the free- 
molecular and transition regimes. High-fidelity CFD 
(LAURA) and low-fidelity Newtonian codes are used 
in hypersonic flow regimes. Finally, both CFD and 
wind tunnel results are employed in the supersonic, 
transonic, and subsonic flow regimes. 

The packaging of the Stardust sample return capsule 
in combination with the light-weight forebody PICA 
heatshield produced a number of c.g. related stability 
issues. These issues were identified with the combined 
trajectory and aerodynamic analyses. For example, 
in the free-molecular and transitional flow regimes, 
a potential for the vehicle to flip around was identi- 
fied. Another issue addressed was the Mach number 
for parachute deployment which is a strong function 
of the transonic and subsonic dynamic stability of the 
sample return capsule. 

Galileo 

The Galileo probe was designed to withstand the 
harshest environment ever encountered by any plan- 
etary probe. 108 One of the major objectives was to 
determine the atmospheric structure of Jupiter. 109 The 
probe entered the Jovian atmosphere at a relative ve- 
locity of 47.5 km/s. and within the first 100 sec. after 
entry, the probe decelerated to less than 1 km/sec. The 
45 deg. blunted-cone probe was designed to withstand 
stagnation point peak-heating rate of 30kW/cm 2 and 
total heat-load of 300kJ/cm 2 . 107 The forebody heat- 
ing was mainly due to shock-layer radiation during 
reentry. These environments were originally defined 
by state-of-the-art computational aerothermodvnamic 
methods of the 1970’s and early 1980’s as represented 
by several papers 160, 161 and reviewed by Howe et. 
ai . 162 

The determination of the structure of the Jovian 
atmosphere from the flight measurements requires ac- 


curate prediction of the the aerodynamic forces in- 
cluding drag coefficient of the probe during entry. 159 
Ballistic range measurement coupled with CFD anal- 
ysis 163 established the drag coefficient in the post- 
and pre-ablative regions. At present the drag co- 
efficient estimation is limited to heurestic modeling 
during the ablation phase. Though the recession mea- 
surement during flight helps in determining the mass 
loss and shape change, we have yet to establish valid 
methodologies in determining the aerodynamic forces 
and moments during the ablative phase. Since accu- 
rate prediction of the aerodynamic forces and moments 
are not only necessary for probe design and trajectory 
determination, these predictions form the basis of “Sci- 
ence Experiments” and it is a challenging task when 
the environment is as harsh as the one encountered by 
the Galileo probe. 

A harder problem during the probe design was the 
heat-shield thickness determination 158 and this poses 
a formidable challenge even today. The ablative heat- 
shield material selected was carbon-phenolic and a 
number of modeling studies were used to determine the 
expected recession of the forebody heat-shield. The 
ablative response modeling of the heat-shield involved 
a coupled thermal environment prediction, including 
radiation and turbulence models with in-depth ma- 
terial response. During the entry, the heat-shield 
experiment allowed the determination of the actual 
surface recessions. The flight measurements compared 
with deign predictions showed the stagnation point re- 
cession to be significantly less, whereas the frustum 
recession far exceeded the predictions. The challenge 
facing our community today is to build and verify In- 
fidelity analysis tools that can accurately predict the 
recession rates for future missions to Jupiter, Europa 
and Neptune. 

STS 

CFD has worked primarily in a post flight, reac- 
tive mode regarding impact on the Shuttle operations 
and flight analysis. Several important studies are cited 
here because they are representative of the role CFD 
plays in this environment. 

An apparent anomaly in the body flap effective- 
ness (relative to pre-flight data book) at high Mach 
number observed in STS-1 has been traced to small 
changes in pressure associated with changes in the 
specific heat integrated over a large, expansion area 
on the windside. 135,164 Thermocouples on early Shut- 
tle Orbiter flights have been used as a critical part of 
the code validation process, 49 including effects of finite 
wall catalysis. The overset grid capability in OVER- 
FLOW has been applied to a transonic, ascent flight 
condition for the complete Shuttle ascent vehicle, in- 
cluding an approximate treatment of real gas effects 
in the plume. 165 A preliminary demonstration of the 
overset methodology for a supersonic flow of a perfect 
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a) Windside 


a) Windside 



b) Leeside 

Fig. 8 B1001 configuration (RLV) 


gas involved a proposed canard at Mach 5.8 on the 
Shuttle. 146 

X-33 and Reusable Launch Vehicle - Phase I 

Numerical simulations of hypersonic flow over pre- 
liminary configurations proposed by Lockheed Mar- 
tin, McDonnell Douglas, and Rockwell for a Reusable 
Launch Vehicle (RLV) and X-33, a technology demon- 
strator for the RLV, were conducted. The process 
for defining the aerothermal characteristics of all three 
configurations by independent CFD teams was essen- 
tially the same. The process, as applied to the Lock- 
heed Martin configuration, is described below. 

The simulations were executed using both chemical 
equilibrium and nonequilibrium gas models. Simula- 
tions were generated over six representative trajectory 
points for descent of the B1001 RLV configuration 
in order to establish traceability of aerothermody- 
namic design issues. Simulations were generated over 
five representative trajectory points for descent of the 
B1001A X-33 configuration. Trajectory points for sim- 
ulation were chosen near peak heating and peak dy- 
namic pressure; other points were selected on the basis 
of convenient anchors for Mach number and angle of 



b) Leeside 


Fig. 9 B1001A configuration (X-33) 

attack variation. Representative surface heating, tem- 
perature and pressure distributions were provided to 
the design team, some examples of which are presented 
here. Procedures for incorporating CFD solutions into 
engineering code (MINIVER) format for subsequent 
use by the thermal design team are also discussed. 

Configurations and Grid 

Configuration with designation B1001, presented in 
Fig. 8, was used in the RLV simulations. It has a ref- 
erence length of 1419.25 inches. Configuration with 
designation B1001A, presented in Fig. 9, was used in 
the X-33 simulations. It has a reference length of 752.2 
inches which is 0.53 of full scale RLV. Both configura- 
tions assume a moment center at 0.66 reference length 
behind the nose. The vehicle geometries are identical 
to scale from the nose to upstream of the wing (hyper- 
vator) root. B1001A is tapered more toward the base 
to reduce base drag as compared to its predecessor. 
It has body flaps on wind and lee sides that termi- 
nate at the cowl trailing edge but extend across most 
of the base lateral dimension. B1001 has no control 
surface preceeding the central base region surrounding 
the aerospike engines. Instead, there is an expansion 
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a) Windside 


a) Windside 



b) Leeside 

Fig. 10 B1001 surface grid (RLV) 


surface in the central region with outboard body flaps 
that extend past the cowl trailing edge. Parts of the 
body flap that extend past the trailing edge of the RLV 
(B1001) were not modeled in Phase I studies because 
body shape had already evolved to the B1001A based 
on wind tunnel test results. A set of simulations at a 
single body flap deflection (50 deg. - 30 deg relative 
to the waterline of the B1001A configuration) will be 
discussed. 

Surface grids for B1001 (RLV) were constructed in 
four sections as shown in Fig. 10. Solutions were gen- 
erated in each section sequentially in a block marching 
mode for B1001. The first section extends from 
the nose to the first terminal plane approximately 10 
inches upstream of the wing (hypervator) root. The 
£ coordinate lines extend from the axis of the nose 
straight back to the terminal plane with roughly equal 
arc length between lines. The coordinate lines extend 
circumferentially around the body. The grid density 
in the first section is (52 x 64) cells. 

The second section continues from the terminal 
plane of the first section to a second terminal plane 
that precedes the vertical tail on the leeside and junc- 
ture of the outboard body flaps and the inboard ex- 



b) Leeside 

Fig. 11 B1001A surface grid (X-33) 


pansion surface. Additional circumferential (/?) coor- 
dinate lines are added in this section to define the wing 
(hypervator) and provide smooth transition of circum- 
ferential arc length in the vicinity of the wing / body 
junction. Some additional coordinate lines were also 
added to provide a smoother transition to the tail and 
bodyflap junction. The section was contructed with 
9 blocks in the circumferential direction and a total 
density of 18 x 116 cells. 

The third section extends from the second terminal 
plane to a third terminal plane at the bodyflap hinge 
line. Additional circumferential density was added to 
define the vertical tail and the outboard bodyflap /in- 
board expansion surface juncture. The section was 
contructed with 12 blocks in the circumferential direc- 
tion and a total density of 18 x 186 cells. 

The final section extends from the third terminal 
plane to the cowl trailing edge. The section was con- 
tructed with 23 blocks in the circumferential direction 
and a total density of 6 x 293 cells. 

Surface grids for B1001A (X-33) were constructed 
in two sections as shown in Fig. 11. Lessons learned in 
the first gridding procedure and a simpler aft geometry 
enabled the simpler grid system. Solutions were gen- 
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Fig. 12 B1001A unstructured surface and volume 

grid (X-33) from FELISA 

erated in each section sequentially in a block inarching 
mode for B1001A. In some cases, solutions were then 
regenerated in a fully coupled mode. 

The first section extends from the nose to the first 
terminal plane approximately 10 inches upstream of 
the wing (hypervator) root. The £ coordinate lines 
extend from the axis of the nose straight back to the 
terminal plane with roughly equal arc length between 
lines. The q coordinate lines extend circumferentially 
around the body. The grid density in the first section 
is (64 x 64) cells. 

The second section continues from the terminal 
plane of the first section to the cowl trailing edge. Ad- 
ditional circumferential (??) coordinate lines are added 
in this section to define the wing (hypervator) and pro- 
vide smooth transition of circumferential arc length in 
the vicinity of the wing / body junction. Some ad- 
ditional coordinate lines were also added to provide a 
smoother transition to the tail and bodyflap junction. 
The section was contructed with 12 blocks in the cir- 
cumferential direction and a total density of 100 x 247 
cells. 

In contrast to the structured grid used by LAURA, 
an unstructured grid used in FELISA analyses is 
shown in Fig. 12. This grid is for an inviscid ap- 
plication, but still reveals some of the requirements 
for relatively fine grid over the nose and in front of 
wing leading edges where the bow shock lies close to 
the body. FELISA has no automatic mechanism for 
collapsing the outer boundary to just in front of the 
bow shock as in LAURA. 

Thermal Analyses 

The most relevant validation data set(s) for X-33 
heating rates include: (1) thermophosphor tests on 
the identical configuration at Mach 10 and Reynolds 
number 500,000 and total temperature 1000K; and (2) 


STS-2 thermocouple data in flight between Mach num- 
bers of 24.3 and 12.86, between Reynolds numbers of 
171,000 to 4,387,000, and between altitudes of 72.4 and 
54.8 km. These tests showed agreement with experi- 
mental data for laminar heating rate over most of the 
windside within 10% with some discrepancies as large 
as 20%. 

The computational aerothermodynamic analyses 
were focused on defining global temperature distribu- 
tions around the RLV and X-33. Thermal analysis 
of the tanks required time dependent data in a read- 
ily accessible format as commonly provided by the 
MINIVER code. 166 The required temporal resolution 
was much finer than the matrix of points considered 
by LAURA alone. The MINIVER code is capable 
of making reasonably accurate estimates of center- 
line heating distributions on vehicles like RLV and 
X-33. However, three-dimensional flow effects ocurring 
off-centerline generally are not well approximated by 
MINIVER analyses without some externally derived 
corrections. The necessary corrections are provided by 
LAURA at off centerline locations at the times defined 
in the CFD matrix. 

Implementation of this procedure in phase I for both 
RLV and X-33 analyses occurred as follows. Heat- 
ing and temperature distributions over the vehicle 
were generated by LAURA and compared with the 
windward centerline results from MINIVER. These 
comparisons established MINIVER as a reasonably 
accurate tool for the geometries and trajectories con- 
sidered here in Phase I studies. Off-centerline values 
are keyed to centerline values of laminar heating rate 
in a relatively dense matrix of computational planes 
through plots of qLaminar / QLaminar, CL as a function 
of spanwise location in the plane. This data is input 
to MINIVER in tabular form. 

Transition to turbulence is assumed to occur for val- 
ues of Res / A/,, between 250 and 300. A transition front 
is defined in LAURA across a computational plane. 
Turbulent heating levels are computed downstream of 
this plane and values of qTurbuient/ QLaminar are de- 
fined using earlier laminar solutions. These turbulent 
to laminar factors are also input into MINIVER in tab- 
ular form. Heating at any point on the body is then 
predicted by MINIVER by computing the windside 
centerline value at the same axial location, multiply- 
ing by an appropriately interpolated value for q/qc.L 
for the spatial location on the body and temporal loca- 
tion along the trajectory, and applying an additional 
correction factor for turbulent flow if the transition 
criteria is exceeded. 

Prediction of the transition criteria by LAURA and 
MINIVER along the windward centerline were in sig- 
nificant disagreement. LAURA predicts the thresh- 
hold transition criteria to occur earlier in the trajec- 
tory than MINIVER. Because these transition criteria 
have historically been derived from engineering code 
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RLV Temperature Contours - 1200 seconds 

Laminar 




a) windside 


RLV Temperature Contours - 1200 seconds 

Laminar Leeside 



Fig. 14 Windside centerline comparisons of 
LAURA and MINIVER temperatures at the 1200 
s trajectory point for the RLV. 

subsequent figures are in degrees Fahrenheit.) Highest 
heating rates occur near the wing root but are only 
slightly higher than the stagnation point heating on 
the nose. 

Solutions at this trajectory point are generated se- 
quentially across four sections of the vehicle. Each 
successive section has more circumferential resolution 
to define the wing, tail, and edges of control surfaces. 
The solution was not run in a fully coupled mode; con- 
sequently, some cosmetic effects of domain interfaces 
are evident in the solution. Reasonably good agree- 
ment with the engineering code MINIVER along the 
windside centerline is demonstrated in Fig. 14. 


RLV Temperature Contours - 1200 seconds 

Laminar 


grid related discontinuities., T °F 



c) side 

Fig. 13 Temperature contours on RLV for fully 
catalytic wall, laminar flow, at 1200 s. 

analyses like MINIVER, Phase I studies proceeded 
using MINIVER criteria. Establishment of a proper 
criteria is a subject of ongoing research. (See related 
discussion in Physical Models.) 

Temperature maps of the vehicle as predicted by 
LAURA for the 1200 s point of the RLV entry are 
presented in Fig. 13. (All temperatures in these and 


Body Flap 

The baseline configuration for B1001A has an ex- 
pansion surface of approximately 20 deg. that leads 
to the cowl trailing edge above the engines and in- 
cludes a stowed body flap. A flap deflection of 50 
deg. (30 deg. into the unexpanded flow) was simu- 
lated to study effects on aerodynamics and heating 
using component isolation of the deflected flap and 
surrounding area. The initial solution came from an 
undeflected flap case. Enhanced grid around the de- 
flected flap was implemented by partitioning the aft 
section into two streamwise sections with 14 blocks 
each. The section immediately preceeding the de- 
flected flap was resolved into 28 streamwise cells by 
149 circumferential cells distributed across 14 circum- 
ferential blocks. The section over the deflected flap was 
resolved into 22 streamwise cells by 173 circumferen- 
tial cells distributed across 14 circumferential blocks. 
Approximately 30 hours of computer time were de- 
voted to obtaining this solution, but the error norm 
never settled down during the simulation. Some ini- 
tial difficulties were attributed to a transient reverse 
flow that set up over the edge of the flap in the exit 
plane. Vacuum boundary conditions 19 were applied to 
survive this transient; nevertheless the solution never 


24 OF 35 


American Institute of Aeronautics and Astronautics Paper 97-2473 





a) windside, fine grid 



b) windside, coarse grid 

Fig. 15 Temperature contours on X-33 for fully 
catalytic wall, laminar flow, at 370 s with a 50 deg. 
bodyflap deflection. 

fully converged. 

Previous experience in computing separated flow in 
front of ramps has indicated that such flows take a long 
time to set up. It appears that the separation bub- 
ble grows slowly with the upstream separation point 
moving further upstream as more mass is entrained. 
Coarsening the grid had been observed to speed up 
the solution process by providing a better initial con- 
dition for a subsequent fine grid solution. In the case 
of X-33, the coarse grid solution (every other point 
deleted) converged very quickly; however, the subse- 
quent fine grid solution (all points restored) would 
rapidly evolve into an apparent unsteady flow within 
the bubble. (The LAURA simulation uses pseudo time 
advancement with large, constant Courant number 
and asynchronous relaxation; consequently, the evo- 
lution of the flow in the bubble is not simulated in 
a time accurate manner.) Embedded vortices appear 
within the bubble but change size and location in the 
simulation. Local hot and cold spots occur in the un- 
steady, fine grid solution (Fig. 15 (a)) that are not 
evident in the steady, coarse grid solution (Fig. 15 (b)). 
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The bodyflap was somewhat hotter in the fine grid 
solution than in the coarse grid though the outboard 
edge of the flap showed equivalent, high radiative equi- 
librium wall temperatures (2400 F). Application of a 
Baldwin-Lomax algebraic turbulence model across the 
separated flow region and the deflected flap also sup- 
pressed unsteady phenomena on the finest grid. 

Simulation of unsteady separation in a hypersonic 
environment on a relatively complex configuration re- 
mains a challenge to state-of-the-art application tools. 
The results (steady vs. unsteady) are sensitive to 
numerical (grid-related) and physical (turbulent vis- 
cosity) dissipation levels. Turbulent viscosity in the 
present test is based on an algebraic model which is 
not appropriate for massively separated flow. Time ac- 
curate simulations involving two-equation turbulence 
models are planned for this case. 

X-33 - Phase II 

In the X-33 Phase I design, the aerothermal en- 
vironment predictions were made for a preliminary 
trajectory (trajectory based TPS design) and config- 
uration. In Phase II, trajectories and configurations 
have evolved to a final design at a pace faster than 
can be accommodated by CFD. Consequently, CFD 
simulations are performed at a number of discrete de- 
sign points described by the variation in Mach number, 
angle-of-attack, and Reynolds numbers. These design 
points are selected to adequately span the design space 
for all possible trajectories. The CFD solution set for a 
given OML (Outer Mold Line) can then be used to de- 
termine the aerothermal heating and heatload for any 
trajectory within the design space quickly through in- 
terpolation with engineering codes. The aerothermal 
environment can be defined in a matter of minutes 
through this approach compared to weeks by the pre- 
vious approach. Once the TPS subsystem is designed 
to a specific trajectory, the TPS margins can then 
be determined quickly by simulating the environment 
for other “off-nominal” trajectories to provide the de- 
sired margin. (A series of papers to be presented at 
the Aerospace Sciences Meeting in January 1998 detail 
the ’’Design-space CFD approach to TPS design” and 
other issues relevant to margins and errors.) 

Forebody to Cowl Trailing Edge 

Both LAURA and GASP codes have been used in 
predicting the aerothermal environment which will be 
used to design the TPS in Phase II. LAURA solu- 
tions are primarily used to define the forebodv, body 
flap and the base region surface temperature at spe- 
cific trajectory points on descent using a high fidelity 
grid. GASP solutions are generated for the forebody 
on a coarser surface grid but on a richer solution point 
matrix, based on complete trajectories and at critical 
design points. In general, for the grids used by the 
two solvers, the GASP and LAURA solutions are in 
good agreement over the majority of acreage on the 
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X33. In fact, the predictions from the two solvers 
have been compared and the difference between the 
two predictions differ no more than 15 F under iden- 
tical freestream conditions using the same grid in the 
nose region. The level of agreement between the two 
solvers is indeed gratifying considering that the numer- 
ical schemes and the procedures are different. Applica- 
tion of two codes by two groups (LAURA at Langley, 
GASP at Ames) at a few overlapping points provides 
an independent, critical check of predictions which en- 
hance confidence in all of the CFD predictions. A brief 
review of the GASP and engineering code HAVOC 167 
application to X33 TPS design follows. More com- 
plete details will be presented in a series of papers to 
be presented in the AIAA Aerospace Sciences Meeting 
in January 1998. 

The use of GASP solver towards the X33 TPS de- 
sign is to generate high-quality, cost-effective solutions 
quickly. A series of grid refinement and grid sequenc- 
ing studies were performed to determine the optimal 
grid points and convergence criteria for a prescribed 
tolerance for surface temperature, which is 25 F or 
less. The total number of grid points required to ac- 
complish this was less than ( 113 * 113 * 65 = 830,000) 
and a single converged reacting-air solution required 35 
hr. on CRAY C-90. A maximum of 10 solutions per 
trajectory were generated to describe the aerothermal 
environment. The integration of the engineering code 
HAVOC coupled with the CFD determined the total 
heat-load and the surface temperature as a function of 
time. 

The aerothermal environment can be defined inde- 
pendent of the trajectory and this approach is cur- 
rently adapted in the X33 program. The design space 
can be discretized in an intelligent manner and a set 
of CFD solutions (approximately 40) can provide the 
aerothermal environment basis set. From these solu- 
tions and using engineering codes such as HAVOC, 
a complete trajectory based solution can be obtained 
in a matter of minutes. Such an approach can be 
integrated with trajectory optimization codes to de- 
termine the optimal trajectory from TPS material per- 
spective. 

Body Flap and Wake 

Wake flow simulations have been implemented for 
X-33 (Phase II) in order to better assess bodyflap effec- 
tiveness and aerothermal loads on the aerospike engine 
on descent. In the present configuration (B1001F), the 
bodyflap extends past the cowl trailing edge into the 
wake. An accurate assessment of bodyflap effective- 
ness must allow for the flow to spill off the sides of 
the flap into the wake. The grids used in Phase II on 
the B1001F are very similar to those generated for the 
B1001A of Phase I. 

Results are presented for a laminar, steady, nonequi- 
librium flow (5 species) on a rj 4 million cell grid 



Fig. 16 Temperature contours in symmetry plane 
of X-33 wake at Mach 10.5 and a = 26 deg. 



Fig. 17 Temperature contours in base region of 
X-33 wake at Mach 10.5 and a = 26 deg. 

distributed across 175 computational blocks. The 
LAURA program requires 228 MW of memory (180 
MW if the Jacobians are stored on disc) and over 
100 C-90 hours for a complete wake flow simulation. 
Extensive use of mesh sequencing and component iso- 
lation procedures were employed in generating these 
solutions. Surface heating levels changed less than 5% 
after the last grid doubling, with most of that change 
concentrated near component edges. 

Figs. 16 and 17 show temperature contours in the 
plane of symmetry and base of the X-33 wake region 
at Mach 10.5 and a = 26 deg. The recompression 
shocks preceeding the wake core are evident in the fig- 
ure. Impingement of the shear layer on the aft end 
of the aerospike is also evident. Surface temperatures 
are evaluated using radiative equilibrium wall bound- 
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ary conditions. The implementation of this boundary 
condition assumed radiation into free space with no 
accommodation given for direct radiation from ad- 
joining surfaces. 168 The engine block temperatures are 
computed with emissivity of 0.12, as compared to 0.8 
for the bodyflap and 0.6 for all other surfaces. Con- 
sequently, the windside engine surfaces respond to a 
given heating rate with higher radiative equilibrium 
wall temperatures than other surfaces. 

Fig. 18 shows the aerodynamic coefficients as a func- 
tion of flap deflection angle measured in the Mach 20 
Helium tunnel at Langley and computed at flight tra- 
jectory points. For zero flap deflection, there is little 
evidence of real gas aerodynamic effects as judged by 
the good comparison between computation at flight 
conditions and ground-based experiment. This trend 
was also observed in Phase I ground based experiments 
when comparing aerodynamics measured in Helium, 
air, and CF 4 . The configuration has a relatively flat 
windside surface as compared to the Shuttle Orbiter 
and shows considerably less sensitivity to real gas ef- 
fects on aerodynamics with undeflected control sur- 
faces. Bodyflap effectiveness as computed for flight 
conditions is stronger than indicated in the experi- 
mental data. This discrepancy is believed to a be 
caused by differences in the separation zone preceeding 
the deflected flap shock strength associated with the 
boundary-layer thickness approaching the flap and ra- 
tio of specific heats. Other relevant validation data 
sets for X-33 aerodynamics are: (1) STS-2 aerody- 
namic data measured in flight between Mach num- 
bers of 24.3 and 12.86, between Reynolds numbers of 
171,000 to 4,387,000, and between altitudes of 72.4 and 
54.8 km. and (2) STS-1 measurement of bodyflap de- 
flection required for trim during a “pitching-moment 
anomaly”. Computed aerodynamic coefficients were 
within experimental uncertainty for the STS-2 data 
points. The computed bodyflap trim angle for STS-1 
was within 10% (1.5 deg.) of the recorded value at a 
= 23. 

Computations that include the wake show flow recir- 
culating back onto the leeside of the vehicle. As with 
the case of separated flow on COMET, 112 extrapola- 
tion outflow boundary conditions on the leeside of a 
vehicle can falsely suggest that the flow stays attached 
to the trailing edge of the vehicle. Misprediction of 
this separation has little impact on hypersonic aero- 
dynamics because of the very low relative pressure on 
the leeside at moderate to high angles of attack. For 
low supersonic to transonic flows, misrepresentation 
of leeside separation can have serious consequences on 
prediction of normal forces and pitching moment. 

There is almost no validation data available for 
the simulation of aerodynamic heating on the engine 
block in the wake. The LAURA code compared well 
to ground based data at Mach 10 for heating levels 
due to free shear layer flow impingement on a sting 





Fig. 18 Aerodynamic coefficients measured in 
ground based experiments and predicted for the 
flight environment as function of bodyflap deflec- 
tion angle for B1001D configuration at a = 36 deg. 
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in the wake of a blunt body using laminar, steady 
flow assumptions. 136 However, transitional effects have 
a strong influence on impingement heating levels in 
the base region and are underpredicted using lami- 
nar flow assumptions. 169 Unlike the bodyflap solution 
from Phase I work, the present wake flow solutions 
showed no indication of unsteady flow. Two-equation 
turbulence models are currently being evaluated for 
application to this problem. 

Leveraging 

Leveraging refers to the use of engineering approx- 
imations, 170 often based on boundary-layer methods, 
to extend (leverage) a limited matrix of CFD solutions 
for better coverage of the parameter space. 

Recent experiences with the Phase I design process 
for X-33 have revealed opportunities for development 
of software to better exploit a limited CFD solution 
matrix. Extraction of only a few parameters at the sur- 
face and at the boundary-layer edge of a CFD solution 
can enable analytic extension of heat transfer solutions 
beyond the baseline matrix. The approach (defined as 
Method 1G 171,172 ) extracts CFD derived quantities at 
the wall and at the boundary layer edge for inclusion 
in a post-processing boundary-layer analysis. It allows 
a designer at a workstation to address two questions, 
given a single CFD solution. (1) How much does the 
heating change for a thermal protection system with 
different catalytic properties than was used in the orig- 
inal CFD solution? (2) How does the heating change 
at the interface of two different TPS materials with an 
abrupt change in catalytic efficiency? The answer to 
the second question is particularly important, because 
abrupt changes from low to high catalytic efficiency 
can lead to localized increase in heating which ex- 
ceeds the usually conservative estimate provided by 
a fully catalytic wall assumption. Design iterations 
are conducted without need of additional CFD runs 
until convergence on a single concept is achieved, at 
which point CFD could be used to provide a final check 
and/or recalibration point. Perhaps more importantly, 
it allows the design team to assess the effect of changes 
in some material properties heating rate without the 
need to rerun archived solutions. 

The front section of the RLV configuration for a 
flowfield simulation at Mach 25, 45 deg. angle of at- 
tack, and 79.6 km altitude is examined. Heating rates 
for both a fully catalytic wall (% = 1 in Eq. 5) and 
a finite-catalytic wall were computed at the respective 
radiative equilibrium wall temperatures. Circumfer- 
ential heating distributions as a function of compu- 
tational coordinate j varying from leeside (j = 1) 
to windside (j = 64) are presented in Fig. 19a for 
the front section. The circumferential cut is from 
the i = 40 plane which lies far downstream from the 
nose and upstream of a wing. The windside centerline 
heating distribution as a function of computational co- 



i 


a) i = 40, circumferential distribution 



b) j = 64, centerline distribution 

Fig. 19 Comparison of CFD heating levels ob- 
tained on front section of RLV with Method 1G 
predictions obtained at 7 = y(T) and 7 = 1 at Mach 
25 and 79.6 km. 

ordinate i varying from the stagnation point ( i = 16) 
to the exit plane of the front section (i = 52) is pre- 
sented in Fig. 19b. Integral boundary-layer corrections 
to CFD baseline heating are generally within 5% of 
computations. Even in the case where significant turn- 
ing of streamlines occurs for flow expanding around 
the side of the vehicle (32 < j < 48), in Fig. 19a the 
integral-boundary-layer extrapolation from the base- 
line CFD computation is an excellent predictor of the 
CFD result at the off-baseline catalysis model. 

Consistency of these predictions over a broader 
range of entry conditions and geometric complexity re- 
mains to be established before integral boundary-layer 
methods (or a related approach) can be used with con- 
fidence in a design mode. 

Boundary-layer codes may extract pressure and ve- 
locity fields from inviscid solutions over complex con- 
figurations to obtain heating distributions. Many more 
solutions per design cycle can be generated with this 
approach (when it is applicable) than with PNS or 
TLNS simulations. The Langley Approximate Three- 
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Fig. 20 Comparison of DPLUR/LATCH and 
LAURA laminar heating in BTU/ft 2 -s on X-33 at 
Mach 11.47 and a = 36.2 deg. 

dimensional Convective Heating (LATCH) 173 code has 
been extensively applied to X-33 and X-34 heat- 
ing analyses using inviscid solutions predominantly 
from DPLUR but also from LAURA and FELISA. 
Interaction with FELISA is complicated at present 
by the requirement to interpolate unstructured grid 
results into the single-block, structured format re- 
quired by LATCH. LATCH combines an axisymmetric 
analog for solution of the general three-dimensional 
boundary-layer equations in a generalized coordinate 
system. Inviscid surface streamlines and pressures are 
interpolated from inviscid solutions; streamline redis- 
tribution is used to obtain uniform coverage of complex 
surfaces. The method developed by Zobv et. al 174 
provides a fast, approximate method for solving the 
boundary-layer equations which has been shown to 
yield accurate results (within ±10% of experimental 
data) for both wind-tunnel and conditions. LATCH 
is currently restricted to single-block, perfect-gas or 
equilibrium-gas applications. Laminar or turbulent so- 
lutions can be generated. While it does not yet handle 
the physical complexity of thermochemical nonequi- 
librium as in codes like the Boundary Layer Integral 
Matrix Procedure (BLIMP) 175 it is applicable to geo- 
metrically complex flows. 

A sample application of LATCH to X-33 heating 
analysis is shown in Fig. 20 and compared to LAURA 
results. Agreement is generally very good. Even 
heating rates on the wing leading edges, where three- 
dimensional effects are exceptionally strong, are within 
20% of the LAURA results. The LATCH solution runs 
in minutes on a workstation. The inviscid DPLUR so- 
lution runs approximately four times faster than the 
viscous LAURA solution. 

Future 

The ongoing industry led design effort for X-33 has 
had a profound effect on the development and applica- 
tion of computational aerothermodynamic tools within 
NASA. The principle drivers in this environment have 


been: (1) the need to respond to evolving configura- 
tions with time consuming gridding requirements; (2) 
the need to produce and evaluate matrices of solutions 
in a tight time frame; and (3) the need to extract, 
interpret, and transfer appropriate subsets of the solu- 
tions in formats that are usable by other members of 
the design team. These drivers have been most valu- 
able in understanding the capabilities and limitations 
of the computational tools that have been developed 
over the past decade. 

Perhaps nothing is more inspiring or motivating 
to a code developer than daily confrontation with 
inadequate software and demanding timelines. In 
some cases, software modifications required to im- 
prove the work environment are simple and incre- 
mental. For example, several utilities (STRIDE, 
ENSEMBLE, BLOX, GRIDSWAP) have evolved in 
LAURA to facilitate grid sequencing, solution se- 
quencing, block marching, component isolation, and 
volume grid restructuring to address X-33 require- 
ments. More flexible interblock boundary conditions 
and routines that apply integral boundary-layer heat- 
ing analyses to CFD solutions were also specifically 
developed in LAURA for X-33. Many other needs 
are also easily identified but not so easily satisfied. 
As noted previously, slow convergence in the wind- 
side boundary layer and in separated regions preceding 
deflected bodyflaps can be partially overcome using 
block marching, component isolation, and physically 
motivated load balancing. However, recent advances 
(predominantly academic in nature) suggest that more 
substantial improvements could be achieved with new 
multigrid methods and pre-conditioning. Incorporat- 
ing these algorithms in a flow environment with strong 
shocks and thermochemical source terms involves re- 
search and verification that are difficult to implement 
in the midst of an intense vehicle design phase. Sim- 
ilar comments can be made regarding grid generation 
issues and unstructured, Navier-Stokes flow solvers. 

Current projects must proceed using the best tools 
on hand. Future projects will require better tools if we 
are to maintain a competitive edge. The challenge is to 
maintain a balance among research and development 
and application activities to serve the present and the 
future. Some thoughts on future needs specifically re- 
lated to computational aerothermodynamics follow. 

Architectures and Algorithms 

CFD underwent a revolutionary change when com- 
puter architectures moved from serial to vector pro- 
cessing. Current computational techniques are based 
on traditional vector processor architectures that have 
not materially changed in the past 20 years. Although 
some effort is under way to port existing codes to a 
parallel environment, there is little work being done 
to exploit the full potential of these systems for hyper- 
sonic CFD. It is generally acknowledged that future 
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increases in computing speed will come through the 
maturation of parallel processing systems. Most cur- 
rent systems are on the (0)100 processors and it is 
not difficult, using current solution algorithms and 
block domain decomposition strategies, to bring the 
full power of one of these systems to bear on a single 
problem. However, is this current approach rational 
or even applicable when these systems are (0)1000 or 
(0)10000 processors. We need to be looking now for 
new and innovative ways to take advantage of these 
systems. 

Convergence Acceleration 

Methods to simplify systems of governing equations 
by identifying partial equilibrium relations via com- 
putational singular perturbation could be explored. 
These ideas have been demonstrated in complex, ther- 
modynamic systems and may have application in more 
conventional fluid mechanics situations. The tech- 
niques fundamentally deal with the question, “Why 
do some parts of the solution or some areas of the 
flow domain converge more quickly than others? Can 
we exploit this behavior to accelerate overall conver- 
gence?” 

Preconditioning applied to the relaxation process 
offers similar potential for improvement, particularly 
in three-dimensional stagnation regions behind strong 
shocks and in the near-wall boundary layer which are 
characterized by near-zero eigenvalues in the Jacobian. 

Application of Multigrid methods for hypersonic 
flows improves convergence times, but appears to fall 
far short of theoretical expectations. We typically see 
a factor of approximately three in speedup at the ex- 
pense of additional memory overhead. Can we get 
large speedups on complex problems with thermo- 
chemical nonequilibrium gas chemistry? Are similar 
approaches possible that are not hindered by presence 
of strong shocks or strong source terms? What savings 
may be associated in applying specialized algorithms 
to strongly elliptic versus strongly hyperbolic or pre- 
dominantly viscous versus predominantly inviscid flow 
domains? What are the best ways to exploit paral- 
lel computing and implement load balancing for these 
applications? Can asynchronous relaxation of various 
flow domains and/or equation sets be exploited for ad- 
ditional convergence acceleration? 

Physical Modeling 

Advances or new approaches in physical modeling 
are sought which offer improved computational effi- 
ciency. These include models for turbulent flow, tran- 
sition from laminar to turbulent flow, energy exchange 
mechanisms, nonequilibrium kinetics, and radiation. 
Can we add sufficient “intelligence” to the algorithm to 
add or delete elements of the physical models on the fly 
where appropriate to reduce computation time while 
maintaining a physically correct simulation? Some as- 
pects of physical modeling introduce new complexities 


in the equation sets. For example, radiative energy 
transport changes our system of PDE’s to an integro- 
differential equation system. The numerical solution of 
these sets generally employs loosely coupled relaxation 
algorithms that require three to five sequential passes 
between equation sets. How can we do it better? 

Inclusion of all significant plasmadynamic effects in 
both weak and strongly ionized flows requires further 
effort. There is evidence 176 that local sound speed and 
possibly drag are affected by low levels of externally 
produced ionization in ways which are not fully under- 
stood or predicted with the current state of physical 
models. 

Surface Definition and Grid Definition 

We need to be able to quickly define configurations. 
We need to be able to quickly alter configurations. 
We need a platform to execute these definitions and 
alterations that is simple (intuitive) to use. Ana- 
lytic definitions (even very complex, patched ones) are 
preferable because they allow for various optimization 
studies. The platform (CAD system) should contain 
a library of shapes, appendages, configurations that 
can be quickly called up and assembled. Operations 
to rescale, twist, bend, stretch elements and auto- 
matically define intersections of elements as they are 
reconfigured should be available. Surface discretiza- 
tion should be carried along as part of this process. 
(If a wing is twisted, the associated surface grid should 
follow the deformation. If a wing is translated along a 
fuselage, the surface grid intersection points on both 
elements should be coordinated. The preservation of 
this information (orientation) is important to facili- 
tate restart solutions.) Structured grid patches or full 
unstructured grid surface disretizations should be al- 
lowed. Automatic filleting of element intersections 
should be accommodated. New elements should be 
input via stereoscopic imaging, discrete data file, or 
analytic definition. Methods for volume grid genera- 
tion/adaption must be advanced in parallel. Volume 
grid generation (structured and unstructured) takes 
too long (hours to weeks) even in cases when volume 
grids over similar configurations already exist. 

In the comments which preceded, it was assumed 
that flow simulation requires discretization of the flow 
domain. Generally, this domain includes the outer 
mold line of the vehicle, a far field boundary, and 
the space in between. Multi-Disciplinary Optimization 
(MDO) analyses including effects of response of the 
vehicle structure to aerodynamic and thermal loads 
(flexing, ablation) are assumed to require additional 
discretization of the structure. At present, there are 
no known methods of analysis that do not require 
discretization. Numerical methods which start with 
an analytic surface definition (likely a very complex, 
patched set of definitions) that generate their own 
surface and volume discretizations on the fly should 
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be investigated. Possibilities for obtaining analytic or 
semianalytic solutions to augment the simulations and 
reduce or eliminate the need for volume discretization 
should be explored (keeping in mind we do not want 
to sacrifice ability to treat complex configurations). 

Gridding strategies in the future will likely exploit 
overset grids and/or a combination of structured and 
unstructured. The state-of-the-art for unstructured 
grid generation on complex configurations is judged 
to be superior to structured grid capability using the 
simple metric of time to generate a complete grid. 
Once one has the complexity of unstructured algo- 
rithm available is there still advantage to maintain- 
ing structured formulation? Preliminary results with 
codes like DPLUR 142 indicate some performance ad- 
vantages with structured grid formulations on SIMD 
machines. If truly multidimensional viscous and invis- 
cid capabilities are present, unstructured approaches 
may be preferable. 
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